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FOREWORD 


This  research  project  of  numerical  techniques  and  solutions  for 
axisymmetrical  viscous  compressible  flows  was  initiated  in  January 
1969  as  a  part  of  a  doctoral  program  in  Mechanical  Engineering  at 
the  University  of  New  Mexico  with  support  by  the  Air  Force.  The 
original  manuscript  was  submitted  to  the  University  of  New  Mexico 
by  Dr.  Kenneth  W.  Smith  as  his  doctoral  dissertation.  Professor 
Victor  J.  Skoglund*  Department  of  Mechanical  Engineeringt  served 
as  the  faculty  advisor  and  collaborated  in  this  research. 

Viscous  compressible  flows  around  blunt  bodies  are  of  primary 
interest  in  a  number  of  relevant  aerospace  applications  including, 
for  example,  missile  design  and  reentry  of  missile  warheads. 
Functional  solutions  for  these  flows  are  inaccurate  because  of 
necessary  mathematical  simplifications  and  detailed  measurements 
of  flow  characteristics  are  diffictilt  and  costly.  In  this  investigation, 
a  complete  set  of  ntunerical  techniques  for  axisymmetrical,  viscous, 
compressible  flows  around  blunt  bodies  is  presented  and  accurate 
solutions  are  obtained. 

Publication  of  this  report  does  not  constitute  Air  Force  approval 
of  the  reported  findings  or  conclusions.  It  is  published  only  for  the 
exchange  and  stimxilation  of  ideas. 


ABSTRACT 


The  purpose  of  this  investigation  was  to  develop  numerical  tech¬ 
niques  for  solving  axisymmetrical,  viscous^  compressible  flow  around 
blunt  bodies.  Solutions  were  limited  to  syscems  of  ideal  gases  with 
laminar  unseparated  bovindary  layers.  The  bow  wave  was  represented 
by  a  moving  discontinviity.  The  remainder  of  the  system  was  repre¬ 
sented  by  second  order  accuracyt  time  dependent  difference  eqxiations. 

In  this  investigation  the  difference  equations  were  derived  from  basic 
time  dependent,  viscous  compressible  flow  equations  that  were  trans¬ 
formed  into  body  related  coordinates.  In  a  development  phase,  many 
numerical  techniques  were  tested  on  digital  computers  before  adopting 
the  ones  that  were  used  in  obtaining  the  results  that  are  presented. 

The  addition  of  stabilizing  terms  to  basic  difference  equations  was  used 
to  achieve  numerical  stability.  Numerical  experiments  v'ere  performed 
to  minimize  the  effect  of  the  stabilizing  terms  on  the  results. 

Solutions  for  a  hemisphere  forebody  were  obtained  at  Mach  2  and 
4  for  inviscid  flow  and  for  several  Reynolds  numbers.  At  Mach  4, 
solutions  were  obtained  for  a  hemisphere-cylinder  for  inviscid  flow 
and  for  a  Reynolds  nximber  of  4000.  Where  possible,  calculated  and 
experimental  results  were  compared.  Their  agreement  was  satisfactory. 
Prior  solutions  of  viscous  compressible  flow  in  the  afterbody  region  of 
blvmt  bodies  were  not  ff.  ond  in  any  publication. 

It  was  concluded  that  accurate  solutions  for  axisymmetrical, 
viscous,  compressible  flows  can  be  obtained  in  the  forebody  and  eifter- 
body  regions  of  blunt  bodies  using  the  time  dependent  numerical  tech¬ 
niques  of  this  investigation.  In  the  forebody  region  of  a  hemisphere- 
cylinder,  approximate  solutions  may  be  obtained  by  solvitig  the  forebody 
system  alone.  Stabilization  of  the  modified  Lax-Wendroff  technique 
seemed  necessary,  but  improvements  in  wave  fitting,  botmdary  and 
digitizing  techniques  might  eliminate  that  need. 
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1.  INTRODUCTION 


1.  1  Motivation 

Analysis  of  the  performance  of  supersonic  vehicles  is  required  in 
their  design.  Their  flow  fields  are  complicated  by  detached  bow  waves 
and  by  complex  flow  fields  behind  the  waves  which  interact  with  viscous 
boundary  layers.  Functional  solutions  for  viscous,  compressible  flow 
about  blunt  bodies  are  inaccurate  because  of  necessary  mathematical 
simplifications.  Accurate  measurements  are  costly  and  difficult, 
particularly  detailed  measurements  of  flow  characteristics.  In  addition, 
most  measurements  are  limited  to  special  systems  and  may  not  apply 
to  a  new  design.  Recently,  numerical  analysis  has  become  a  valuable 
supplement  to  other  types  of  analyses  and  experiments.  When  better 
numerical  techniques  are  developed  and  larger  and  faster  computers 
become  available,  they  may  become  the  main  analytical  tools  of  designers 
of  supersonic  vehicles.  Numerical  techniques  have  been  developed  and 
solutions  reported  for  viscous,  compressible,  supersonic  flow  in  the 
forebody  region  of  a  two  dimensional  body.  In  this  investigation,  a 
complete  set  of  numerical  techniques  was  developed  for  axisymmetrical, 
viscous,  compressible  flow  around  blunt  bodies. 

1.2  Purpose 

The  purpose  was  to  develop  numerical  techniques  for  axisymmetrical 
viscous,  compressible  flows  in  both  the  forebody  and  afterbody  region  of 
blunt  bodies. 
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1.3  Scope 

The  investigation  included  development  of  numerical  techniques 
and  solutions  of  axisymmetrical,  viscous,  compressible  flows  around 
hemispheres  and  hemisphere-cylinders.  The  solutions  were  restricted 
to  systems  of  ideal  gases  with  laminar  unseparated  boundary  layers. 

In  the  development  phase,  the  time  dependent  differencing  technique 
of  Lax  and  Wendroff  [45,  1960]^  was  used  to  represent  the  basic  flow 
equations  which  were  expressed  in  terms  of  body  related  coordinates. 

The  detached  bow  wave  was  represented  by  a  moving  discontinuity.  It 
was  coupled  to  the  digitized  field  by  a  wave  fitting  technique.  The  addi¬ 
tion  of  stabilizing  terms  to  the  difference  eouations  was  used  to  achieve 
numerical  stability.  Numerical  experiments  were  used  to  minimize  the 
effect  of  stabilizing  terms  'n  results. 

Solutions  were  obtained  for  hemispheres  at  Mach  2  and  4  for 

3  5 

inviscid  flows  and  Reynolds  numbers  from  10  to  10  .  Inviscid  flow- 
results  agreed  with  measurements  of  Baer  [4,  1961].  Computed  veloc¬ 
ities  in  the  boundary  layer  agreed  with  measurements  oi  Wells  and 
Blumer  [75,  19683„  Solutions  were  also  obtained  for  a  hemisphere- 
cylinder  for  Mach  4  for  inviscid  flow  and  a  Reynolds  number  of  4000. 

A  steady  solution  for  a  hemisphere-cylinder  was  used  as  the  initial 
condition  for  a  subsequent  solution  in  which  the  stabilizing  terms  were 
zero. 

^  Numbers  in  brackets  []  designate  references  listed  at  the  enC  of 
this  report. 
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The  "Review"  section  describes  previous  research  that  was  pertinent 
to  this  investigation.  The  "Theory"  section  presents  the  basis  of  the 
adopted  numerical  techniques.  The  section  on  "Numerical  Techniques" 
presents  the  equations  that  were  the  basis  of  computer  programs.  The 
"Computations  and  Results"  section  describes  computer  programs  and 
the  data  and  results  of  this  investigation. 


2.  REVIEW 


2.  1  Introduction 

The  procedure  used  in  the  "Review"  was  to  first  survey  titles  and 
secure  abstracts  of  literature  pertinent  to  solutions  of  blunt  body  systems. 
From  these  titles  and  abstracts,  literature  was  selected  for  further 
study.  Using  the  literature  that  was  collected,  different  methods  of 
solving  for  the  flow  about  blunt  body  systems  were  studied.  From  this 
study,  the  explicit  time  dependent  method  was  selected  as  the  method 
to  be  used  in  this  investigation.  The  purpose  of  this  "Review"  section 
is  to  summarize  Information  that  was  useful  in  this  investigation.  Topics 
that  are  covered  are  methods  for  obtaining  Initial  values,  exph'i',  time 
dependent  differencing  techniques,  wave  fitting,  prior  solutions  and 
experimental  results  for  blunt  body  systems. 

Studies  of  the  literature  collected  were  focused  on  the  advantages 
and  disadvantages  of  available  methods  for  solving  axisymmetrical, 
viscous,  compressible  flows.  Without  high  speed  computers,  the  extent 
of  early  computations  was  limited.  One  of  the  common  methods  used  to 
solve  for  viscous  flow  around  a  blunt  body  system  is  to  divide  the  shock 
layer  into  an  inviscid  and  viscous  region.  Basic  textbooks  on  boundary 
layer  theory  by  Schlichting  [63,  1966]  and  Rosenhead  [42,  1963]  deal 
with  the  Prandtl  boundary  layer  theory  in  which  the  inviscid  solution  is 
matched  with  the  viscous  flow  next  to  solid  surfaces.  In  solving  axisym- 
metrical,  viscous,  compressible  flows  early  methods  required  many 
assumptions  and  results  were  only  approximate.  Recent  availability 


of  high  speed  digital  computers  has  permitted  more  accurate  solutions 
in  which  functional  and  numerical  methods  are  combined.  The  literature 
on  these  methods  is  extensive.  A  few  examples  are  given  in  [32,  1955; 
41,  1959;  23,  1964;  28,  1968;  22,  1968;  and  25,  1969]  in  which  an 
inviscid  flow  solution  provides  boundary  conditions  for  the  boundary 
layers.  Functional  methods  Involve  assumptions  and  provide  only 
approximate  solutions  of  the  inviscid  flow.  Hayes  and  Probstein  [38, 
1966]  deal  with  methods  of  solving  inviscid  flow  systems. 

More  accurate  solutions  of  viscous  compressible  flows  are 
obtained  by  using  finite  differences  to  approximate  the  governing  partial 
differential  equations  of  a  system.  In  addition  to  the  solution  accuracy 
gained  by  using  numerical  methods,  flexibility  is  improved  in  repre¬ 
senting  boundary  conditions.  Numerical  techniques  were  selected  for 
use  in  this  investigation  because  of  their  accuracy  and  flexibility. 
Numerical  techniques  and  solutions  are  described  in  sections  2.  3  and 
2.5. 

In  the  study  of  numerical  techniques,  attention  was  focused  on 
their  advantages  and  disadvantages  for  blunt  body  systems  including 
differencing  and  initial  and  boundary  conditions.  After  studying  the 
literature,  a  class  of  explicit,  time  dependent,  finite  difference  tech- 
x.iques  was  selected  to  be  used  in  this  Investigation.  The  reason  for 
this  selection  was  that  implicit  methods  involve  simultaneous  solutions 
of  all  nodal  data  points  and  are  more  difficult  to  program  for  a  digital 
computer  than  the  straightforward  explicit  finite  difference  approach. 
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For  some  systems,  implicit  methods  require  less  computer  time,  but 
implicit  methods  have  not  been  successfully  applied  to  viscous  compress- 
ible  flows.  For  viscous  flow  around  a  blunt  body,  coordinate  transforma¬ 
tions  are  often  used  which  complicate  the  governing  partial  differential 
equation.  The  explicit  finite  difference  technique  lends  Itself  to  simple 
treatment  of  these  complicated  equations.  Therefore,  an  explicit,  time 
dependent,  differencing  technique  was  selected  because  of  its  greater 
flexibility  and  mathematical  simplicity.  Another  advantage  of  finite 
differencing  partial  differential  equations  in  time  dependent  form  is 
that  these  equations  are  hyperbolic  regardless  of  Mach  number.  This 
is  especially  advantageous  for  viscous  compressible  flows  where  there 
are  regions  of  mixed  subsonic  and  supersonic  flow. 

In  the  explicit  time  dependent  technique  selected  for  use  in  this 
investigation,  dependent  functionals,  such  as  velocity  and  temperature, 
are  determined  only  at  a  finite  number  of  locations  which  are  called 
nodes.  A  solution  is  started  with  initial  values  specified  at  all  nodes. 
Later  values  are  calculated  by  repeated  application  of  the  difference 
equations  that  are  analogs  of  the  governing  partial  differential  equations 
for  specific  time  increments.  Boundary  conditions  must  also  be  repre¬ 
sented  in  finite  difference  form.  Shock  waves  in  the  system  may  be 
represented  by  the  difference  ex^uations  which  spread  the  wave  over 
several  nodes.  Waves  may  also  be  represented  as  discontinuities  in 
which  the  Rankine-Hugoniot  equations  apply. 

Experimental  results  of  previous  investigations  were  studied  for 
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comparison  with  those  of  this  investigation.  Only  steady  flow  results 
from  wind  tunnel  experiments  were  available.  A  special  effort  was 
made  to  obtain  results  for  2  s  M  s  4  where  gas  properties  are  related 
by  ideal  gas  equations. 

2.2  Initial  Values 

2.2.  1  Initial  Bow  Wave.  Bow  wave  coordinates  are  required  to  specify 
initial  values.  Moretti,  et  al.  [51,  1968;  52,  1968;  53,  1968]  used 
relatively  crude  methods  to  initially  estimate  bow  wave  coordinates  for 
an  Inviscid  solution.  Viscous  systems  require  small  nodal  and  time 
incren.ents  which  result  in  long  computing  times  for  convergence. 
Therefore,  a  more  accurate  determination  of  the  initial  bow  wave 
coordinates  was  desirable  in  the  subject  investigation.  An  accurate 
method  of  predicting  bow  wave  coordinates  was  given  by  Moeckel  [50, 
1949]  and  Love  [16,  1957].  Their  techniques  involve  a  combination  of 
empirical  and  functional  methods  which  are  relatively  simple  to  code 
for  a  digital  computer. 

2.2.2  Boundary  Layer.  A  method  for  calculating  laminar  compressible 
boundary  layers  of  a  blunt  body  was  needed  in  determining  initial  values. 
Many  methods  were  reviewed  and  rejected  as  either  being  too  compli¬ 
cated  or  too  restrictive  in  flow  conditions  [32,  1955;  23,  1964;  27,  1953; 
14,  1949].  The  Cohen  and  Reshotko  [19,  1955]  results  compared  favor¬ 
ably  with  the  experimental  boundary  layer  velocities  of  Wells  and  Blumer 
[75,  1968].  Cohen  and  Reshotko  present  tables  of  a  two  dimensional 
solution  for  Prandtl  Number  =  1  and  constant  body  temperature. 
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Pressure  gradients  are  allowed  from  the  infinitely  favorable  to  the 
adverse  gradients  of  separation.  Results  are  given  for  body  surface 
temperatures  from  absolute  zero  to  twice  the  free  stream  stagnation 
value.  Recently,  their  tables  have  been  extended  by  Christian,  Hankey 
and  Petty  [15,  1970].  The  tables  of  both  references  involve  the  assump 
tlon  that  the  Inviscid  velocity  at  the  edge  of  the  boundary  layer  is 

u^  =  rx""  (1) 

Here  F  and  m  are  constants  and  X  is  the  distance  along  the  body  sur¬ 
face.  In  utilizing  the  tables,  segments  of  the  boundary  layer  were 
matched  to  equation  1. 

2.  3  Differencing  Techniques 

2.3. 1  Introduction.  Basic  theorems  for  stability  and  convergence  of 
difference  equations  were  derived  by  Courant,  Fredricks  and  Lewy  [20, 
1928].  Von  Neumann  [73,  1944]  suggested  that  hyperbolic  partial 
differential  equations  representing  mixed  subsonic  and  supersonic 
flows  could  be  solved  by  finite  difference  equations.  Von  Neumann 
and  Richtmyer  [74,  1950]  added  a  stabilizing  term  to  their  difference 
equations  to  stabilize  the  numerical  solution.  Stabilizing  terms  have 
been  developed  which  do  not  seriously  degrade  accuracy  except  in  the 
vicinity  of  shock  waves.  Lax  [44,  1954]  systematically  developed 
differencing  techniques  applicable  to  strong  shock  waves,  but  the 
waves  are  spread  over  several  nodes.  Rusanov  [61,  1962]  minimized 
the  stabilizing  term  used  by  Lax  to  improve  accuracy.  However, 
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Emery  [30,  1967]  reported  that  computations  using  Rusanov's  technique 
were  unstable  after  long  time  periods.  Aungier  [1,  1970]  developed 
simple  stabilizing  terms  as  a  variation  of  Lax's  method  and  his  results 
for  inviscid  flows  compared  very  favorably  with  those  of  experiments 
and  with  those  of  the  method  of  characteristics.  Lax  and  Wendroff  [45, 
I960]  reported  a  differencing  technique  of  second  order  accuracy  which 
has  been  widely  used  both  in  its  original  and  modified  forms.  Other 
time  dependent  differencing  techniques  have  been  introduced  by  Godunov 
[35,  1959],  MacCormack  [47,  1969]  and  Crocco  [21,  1965]. 

After  studying  the  available  differencing  techniques,  the  one  of 
Lax-Wendroff  was  selected  for  differencing  the  field  equations  in  this 
investigation.  It  was  selected  because  successful  applications  have  been 
reported  for  both  inviscid  and  viscous  flow  systems. 

2.  3.2  Lax-Wendroff  Differencing  Technique.  Lax  and  Wendroff  trun¬ 
cated  a  Taylor  expansion  to  yield  difference  equations  of  second  order 
accuracy.  The  technique  is  described  in  detail  in  section  4.  3.2. 

Burstein  [12,  1965]  used  the  Lax-Wendroff  technique  to  solve 
for  the  inviscid,  two  dimensional  flow  over  a  body  with  a  square  nose. 

In  his  solution,  instabilities  occurred  near  the  detached  shock  and 
sonic  line.  Burstein  eliminated  the  instabilities  by  adding  stabilizing 
terms  to  the  Lax-Wendroff  difference  equations.  Lapidus  [43,  1967] 
calculated  the  characteristics  of  Inviscid  flow  over  a  two  dimensional 
cylinder  using  a  modified  version  of  the  Lax-Wendroff  technique. 

Lapidus  also  experienced  difficulties  in  stabilizing  his  solution. 
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However,  he  developed  stabilizing  terms  which  were  less  complicated 
than  those  of  Bur  stein,  and  he  succeeded  in  obtaining  a  solution.  The 
stabilizing  technique  of  Lapidus  is  described  in  section  4.  5.2. 

Skoglund,  Cole  and  Staiano  [65,  1967]  developed  methods  for 
solving  for  the  interaction  of  an  oblique  shock  wave  with  a  laminar 
boundary  layer.  Additional  stabilizing  terms  were  not  necessary  to 
achieve  numerical  stability  for  that  system  using  the  Lax-Wendroff 
technique.  Later,  Skoglund  and  Gay  [66,  1969]  extended  the  work  to 
include  separation  of  the  boundary  layer.  Using  the  techniques  of 
Skoglund,  et  al.  [65,  1967],  instabilities  occurred  in  the  separated 
region  near  the  edge  of  the  boundary  layer.  Numerical  stabilization 
was  accomplished  by  adding  stabilizing  terms  that  were  derived  from 
those  of  Lapidus  [43,  1967]. 

Richtmyer  and  Morton  [59,  1967],  as  well  as  Lapidus  [43,  1967], 
proposed  a  two-step,  Lax-Wendroff  technique.  Values  are  obtained  at 
time  t  using  the  first  order  technique  of  Lax  [44,  1954].  Values 

at  t  =  t+At  are  calculated  using  centered  time  differences  based  on 
values  at  t  Erdos  and  Zakkay  [31,  1969]  used  the  two-step, 

Lax-Wendroff  technique  for  solving  an  inviscld  two  dimensional  flow- 
in  the  near  wake  region.  They  added  a  stabilizing  term  to  the  differ¬ 
ence  equations  of  their  second  step  in  order  to  obtain  a  stable  solution. 
2.3.3  Aungler  Differencing  Technique.  Aungier  [1,  1970;  2,  1971; 

3,  1968]  developed  a  version  of  the  Lax  technique  to  solve  for  inviscld 
compressible  flow  about  blunt  bodies.  His  differencing  technique  is 
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described  in  section  4.  3.  2.  The  stabilizing  term  developed  by  Aungicr 
is  described  in  section  4.  5.  3. 

2.4  Wave  Fitting 

As  indicated  by  references  [30,  1967;  39,  1954;  45,  1967;  65,  1967; 
66,  1969;  74,  1950],  representation  of  shock  waves  by  difference  equa¬ 
tions  results  in  spreading  the  wave  over  several  nodes.  For  the  blunt 
body  system,  it  seems  better  to  represent  the  bow  wave  as  a  disconti¬ 
nuity  which  satisfies  the  Rankine-Hugoniot  equations.  The  latter  tech¬ 
nique  is  called  wave  fitting.  An  important  problem  of  this  technique  is 
the  coupling  of  the  wave  to  the  digiti..ed  field  which  is  represented  by 
difference  equations.  This  technique  was  described  by  Richtmyer  [64, 
1961].  Moretti  [52,  1966;  53,  1968;  and  51,  1968]  used  a  variation  of 
Richtmyer’s  technique  for  blunt  body  systems.  His  method  is  explained 
in  detail  in  [53,  1968]. 

2.  5  Numerical  Solutions  of  Blunt  Body  Systems 

Solutions  of  blunt  body  systems  using  time  dependent  numerical 
techniques  became  feasible  with  the  advent  of  high  speed  digital  com¬ 
puters-.  An  early  paper  of  Bursteln  [12,  1965]  was  on  inviscid  compress¬ 
ible  flow  over  a  two  dimensional,  flat-nosed,  blunt  body.  Free  stream 
values  were  used  for  initial  conditions,  and  the  initial  bow  wave  was  at 
the  body  surface.  The  bow  wave  was  represented  by  difference  equa¬ 
tions  and  was  spread  over  several  nodes.  Burstein  did  not  describe 
the  techniques  that  he  used  for  the  body  surface  and  at  the  downstream 
boundary.  However,  in  an  earlier  paper  [13,  1964],  Burstein  used  a 


reflection  technique  at  the  body  surface  and  a  backward  differencing 
technique  at  the  downstream  boundary.  Reflection  techniques  involve 
strings  of  nodes  along  the  body  surfaces  and  within  the  solid  body. 
Absolute  values  at  the  nodes  inside  the  body  are  set  equal  to  the  values 
at  the  first  string  of  nodes  outside  of  the  body  surface.  Burstein  used 
a  variation  of  the  two-step,  Lax-Wendroff  differencing  technique  for 
points  in  the  field.  Numerical  instabilities  occurred  in  the  stagnation 
region,  near  the  bow  wave  and  near  the  downstream  boundary.  The 
solution  was  stabilized  by  adding  stabilizing  terms  to  the  differencing 
equations. 

Bohachevsky  and  Rubin  [8,  1966]  used  a  Lax  differencing  tech¬ 
nique  [44,  1954]  to  solve  for  the  nonequilibrium  invlscld  flow  over  a 
variety  of  two  dimensional  and  axisymmetrical  bodies.  The  grid  system 
was  extended  well  beyond  the  expected  bow  wave  into  the  free  stream. 

The  outer  and  downstream  boundaries  were  treated  in  the  following 
different  ways:  (1)  first  derivatives  zero,  (2)  second  derivatives  zero, 

(3)  reflection  technique  and  (4)  the  boundary  was  left  free.  The  results 
were  approximately  the  same  for  all  of  these  techniques.  They  concluded 
that  the  choice  was  unimportant  as  long  as  the  flow  at  the  boundary  was 
supersonic.  Body  surface  nodes  were  treated  using  the  reflection  tech¬ 
niques  of  Burstein  [13,  1964].  Shockwaves  were  represented  by 
difference  equations.  Bohachevsky  and  Rubin  did  not  report  any 
trouble  with  numerical  instabilities.  However,  their  results  were 
inaccurate  in  the  stagnation  region.  Bohachevsky  and  Rubin  ascribed 
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the  errors  to  difficulties  in  differencinf  in  spherical  coordinates. 

Later,  Bohachevsky  and  Mates  [7,  1966]  extended  the  solution 
to  include  angles  of  attack.  Additional  problems  were  increased  com¬ 
puter  storage,  further  decrease  in  accuracy  and  large  computer  outputs. 
DeJarnette  [24,  1966]  also  used  the  Lax  technique  to  calculate  invisc*  ' 
nonequilibrium  flow  over  a  blunt  body.  DeJarnette  limited  his  solut. 
to  supersonic  regions  and  relied  on  some  other  solution  to  provide 
initial  values  downstream  of  the  sonic  line.  Using  a  modified  forward 
differencing  technique  at  the  body  surface,  he  reported  that  results  were 
more  accurate  than  those  of  the  reflection  technique.  However, 
DeJarnette 's  results  show  large  perturbations  downstream  of  the  sonic 
line. 

Although  not  reported  by  Bohachevsky  and  Rubin  [8,  1966], 
according  to  an  analysis  by  Moretti  and  Abbett  [52,  1966],  the  required 
computer  running  time  to  calculate  the  flow  around  a  step  was  approxi¬ 
mately  four  hours  on  an  IBM  7094.  The  number  of  nodes  required  to 
provide  good  resolution  of  the  bow  shock  was  3588  and  results  after 
650  cycles  were  presented  in  [8,  1966];  however,  Moretti  and  Abbett 
maintained  that  in  their  opinion,  convergence  still  had  not  occurred. 

Moretti  and  Abbett  represented  the  bow  wave  by  a  discontinuity. 
Since  resolution  near  the  wave  was  not  a  problem,  the  number  of  nodes 
in  the  field  was  reduced,  and  with  fewer  nodes,  computational  times 
were  greatly  reduced,  Moretti  and  Abbett  calculated  inviscid  flows 
over  the  forebody  of  various  shapes.  Initial  values  were  obtained  by 
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assuming  a  parabolic  bow  wave  shape  with  the  standoff  distance  being 
calculated  on  the  basis  of  curvature.  Initial  values  at  the  stagnation 
point  were  calculated  by  assuming  isentropic  flow  behind  the  bow  wave. 
Interpolations  were  used  for  values  in  the  region  between  body  surface 
and  bow  wave.  Linear  extrapolation  was  used  at  the  downstream 
boundary.  Values  at  the  body  surface  were  calculated  using  the  method 
of  characteristics.  Mapping  of  the  field  into  a  rectangle  was  accom¬ 
plished  using  a  simple  coordinate  transformation.  Moretti  and  Abbett 
used  a  Lax-Wendroff  differencing  technique.  In  order  to  speed  conver- 

gZf 

gence,  the  — r  term,  was  arbitrarily  multiplied  by  2.  Accurate  results 

at"^ 

were  obtained  in  the  forebody  region  for  two  dimensional  and  axlsym- 
metrical  flows.  Computer  times  varied  from  15  seconds  to  6  minutes 
on  an  IBM  7094  computer,  depending  on  resolution.  This  vast  improve¬ 
ment  over  the  computational  times  of  Bohachevsky  and  Rubin  is  due  to 
a  reduction  in  the  number  of  nodes. 

The  foregoing  report  has  served  as  a  foundation  for  a  series  of 
reports  by  Moretti  and  others.  Moretti  and  Blelch  [53,  1968]  extended 
the  earlier  techniques  and  obtained  solutions  for  three  dimensional 
inviscld  flow  around  a  blunt  body.  Their  techniques  were  identical  to 
those  employed  by  Moretti  and  Abbett  except  that  the  need  for  a  con¬ 
vergence  term  was  not  mentioned.  A  typical  running  time  was  thirty 
minutes  on  an  IBM  7094  computer.  A  maximum  of  594  nodes  was  used. 
From  the  results  presented,  it  appears  that  convergence  occurred 
within  300  to  500  cycles.  Results  were  limited  to  the  forebody  region 
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because  of  limited  computer  storage.  Later,  Moretti  [51,  19o8] 
improved  convergence  and  accuracy  of  the  results  of  [53,  1968]  by 
using  spherical  coordinates.  This  dependence  of  stability,  convergence 
and  accuracy  on  the  coordinate  system  confirmed  the  findings  of 
Bohachevsky  and  Rubin. 

Moretti  and  Salas  [54,  1969]  extended  the  technique  of  Moretti 
and  Abbett  [52,  1966]  to  include  viscosity  and  thermal  conductivity. 

A  nonlinear  coordinate  transformation  was  used  to  increase  resolution 
of  the  boundary  layer.  The  calculations  of  Moretti  and  Salas  were 
limited  to  the  forebody  region  of  a  two  dimensional,  circular  cylinder. 
They  used  Lax-Wendroff  differencing  techniques  in  polar  coordinates. 

In  all  cases,  body  temperature  and  viscosity  were  constant  and  Mach 
number  was  equal  to  4.  Reynolds  numbers,  that  were  referred  to  free 

2  5 

stream  conditions  and  body  radius,  ranged  from  10  to  10  .  A  typical 

computing  time  v,as  4  minutes  for  1000  cycles  on  a  CDC  6600  computer. 

The  number  of  cycles  varied  from  560  for  Re  ^  5000  to  1500  for 

Re  >  5000.  In  that  time,  convergence  was  not  yet  complete,  but 

Moretti  and  Salas  considered  the  results  acceptable.  Later,  for  the 

same  system,  Moretti  and  Salas  [55,  19/0]  abandoned  the  Lax-Wendroff 

differencing  technique  in  favor  of  the  predictor-corrector  one  of 

MacCormack  [47,  1969].  An  equation  written  in  terms  of  ^  was  used 

ot 

instead  of  the  continuity  equation.  Even  though  coding  was  simplified, 
the  reduction  in  computing  time  was  not  significant.  Although  no  insta¬ 
bilities  were  reported,  Moretti  stated  in  a  private  communication  that 
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instabilities  occurred  if  the  flow  u/as  not  accelerating  at  the  down¬ 
stream  boundary. 

Scala  and  Gordon  [62,  1968]  obtained  a  solution  for  viscous 
compressible  flow  around  a  two  dimensional  cylinder  usln«»  an  explicit 
time  dependent  technique.  Their  technique  involved  a  slight  variation 
of  the  Crocco  differencing  technique  [21,  1965],  The  bow  wave  was 
represented  by  finite  difference  equations.  The  total  number  of  com¬ 
putational  nodes  was  large  in  order  to  achieve  adequate  resolution  in 
the  boundary  layer.  In  one  solution  case,  637,000  nodal  computations 
required  twenty  hours  of  IBM  7094  computer  time. 

Godunov  [35,  1959]  developed  a  time  dependent  technique  which 
is  supposedly  an  optimum  combination  of  characteristic  and  difference 
equations.  Godunov  represented  tne  bow  wave  by  a  discontinuity;  how¬ 
ever,  in  his  calculations  of  the  inviscid  flow  over  a  blunt  body  [36,  1961], 
the  results  were  inaccurate  near  the  stagnation  point  and  the  bow  wave. 
Masson,  Taylor  and  Foster  [49,  1969]  deduced  that  Godunov's  treatment 
of  the  bow  wave  and  body  surface  was  Incortect.  They  obtained  much 
better  results  using  a  slight  variation  of  Moretti's  technique  for  the  bow- 
wave  and  body  surface.  However,  their  results  at  and  near  the  stagna¬ 
tion  point  were  still  in  error.  Aungier  [1,  1970]  also  concluded  that 
Godunov's  technique  does  not  properly  represent  physical  conditions 
along  stagnation  streamlines  of  axisymmetric  flows. 

Lapidus  [43,  1967]  solved  for  the  characteristics  of  inviscid 
flow  over  a  two  dimensional  cylinder.  Initial  values  including  bow  wave 
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coordinates  were  taken  from  the  results  of  Swenson  [69,  1964],  An 
unusual  coordinate  transformation  was  used  to  obtain  a  rectangular 
field.  The  bow  wave  was  represented  by  difference  equations.  A  two- 
step,  Lax-Wendroff  differencing  technique  was  used  for  interior  nodes. 
Instabilities  occurred  when  linear  extrapolation  was  used  at  the  down¬ 
stream  boundary  in  early  stages  of  the  computation.  Computations  were 
stabilized  by  using  an  arbitrary  technique  for  the  initial  500  cycles. 

After  500  cycles,  linear  extrapolation  was  satisfactory.  Instabilkies 
also  occurred  near  the  stagnation  point  and  near  the  bow  wave  that 
were  similar  to  these  reported  by  Burstein  [12,  1965],  Lapidus  used 
a  simpler  stabilizing  term  than  the  one  used  by  Burstein.  The  difference 
of  his  results  is  as  much  as  30%  from  the  more  accurate  ones  of  Swenson, 
even  though  Swenson's  technique  was  used  to  calculate  initial  values. 

Aungier  [1,  1970]  assumed  that  the  initial  bow  wave  was  very 
close  to  the  body  and  that  its  shape  was  the  same  as  that  of  the  body. 
Rankine-Hugoniot  relationships  and  isentropic  flow  equations  were  used 
to  obtain  initial  values  in  the  field.  Linear  extrapolation  was  used  at 
the  downstream  boundary.  Forward  differencing  in  terms  of  body 
related  coordinates  was  used  at  the  body  surface.  The  field  was  seg¬ 
mented  and  a  steady  state  was  achieved  in  one  segment  before  computing 
the  next  segment.  The  bow  wave  was  represented  by  a  discontinuity. 
Aungier  was  able  to  carry  his  time  dependent  inviscid  solutions  into  the 
afterbody  region.  He  reported  that,  within  his  knowledge,  his  time 
dependent  solution  was  the  first  successful  one  for  the  afterbody  region. 
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After  deciding  to  use  explicit,  time  dependent  differencing  tech¬ 
niques  in  this  investigation,  other  numerical  methods  for  blunt  body- 
systems  were  of  little  Interest  and  are  described  only  briefly.  Van 
Dyke  [72,  1958]  gave  a  good  survey  of  the  early  indirect  Inviscid 
methods  in  which  the  bow  wave  shape  was  assumed  known  and  body 
shape  was  then  part  of  the  solution.  A  more  accurate,  but  complicated, 
solution  of  a  two  dimensional,  inviscid,  blunt  body  system  was  given  by 
Swenson  [69,  1965].  Inouge  and  Lomax  [39,  1962]  solved  for  the 
inviscid  flow  over  several  blunt  body  shapes  using  an  Indirect  method 
in  the  forebody  region  where  the  Mach  number  was  less  than  1.03. 

This  solution  was  then  used  as  a  starting  point  for  the  method  of  char¬ 
acteristics  for  the  remainder  of  the  field.  Calculated  body  pressures 
agreed  closely  with  experiment  in  the  forebody  and  afterbody  regions. 
Slight  deviations  of  calculated  and  experimental  results  occurred  in  the 
vicinity  of  the  forebody-afterbody  junction, 

Dorodnitsyn  [26,  1957]  proposed  a  numerical  method  for  non¬ 
linear  flow  equations  which  he  called  the  method  of  integral  relations. 
This  method  has  been  used  by  many  investigators  for  inviscid  flow 
[71,  I960;  70,  1963;  11,  1964;  6,  1965].  The  method  has  an  advantage 
of  being  direct,  so  that  the  bow  wave  coordinates  are  a  part  of  the  solu¬ 
tion.  Unfortunately,  the  method  becomes  extremely  complex  as  the 
resolution  in  the  shock  layer  is  increased  and,  therefore,  seems  imprac¬ 
tical  for  viscous  systems. 
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2.  6  Experimental  Results 

A  large  amount  of  experimental  data  has  been  collected  for  the 
flow  over  blunt  bodies.  However,  reports  containing  data  of  both  the 
body  surface  pressures  and  coordinates  of  the  bow  wave  are  scarce. 
Baer  [4,  1961]  conducted  wind  tunnel  experiments  on  an  AGARD  Model 
E  hemisphere  cylinder  configuration.  Free  stream  Mach  numbers 


ranged  from  2  through  S  while  Reynolds  numbers  varied  from  .  17x10 
.6 


to  .51x10  per  inch.  In  this  report,  coordinates  of  the  bow  wave  are 
specified.  In  many  other  reports,  only  unsealed  schlieren  photographs 
are  available. 

Inouge  and  Lomax  [39,  1962]  summarized  experimental  results 
of  Kendall  [40,  1959],  Kubota  [42,  1957]  and  Baer  [4,  1961].  Graphs 
of  the  bow  wave,  as  well  as  body  surface  pressure  distributions,  for 
a  hemisphere  at  M  =  4.  76,  hemisphere  cylinder  at  M  =  7.  7,  sphere 
cone  at  M  =  4.  95,  and  a  blunt  ellipsoid-cylinder  at  M  =  5. 12  are  pre¬ 
sented. 


Pressure  distributions  over  hemisphere-cylinder  bodies  were 


given  by  Reichle  [56,  1962]  for  Mach  numbers  varying  from  .4  to  5.0 


and  Reynolds  numbers  from  ,36x10^  to  1.2x10^.  Unfortunately,  bow 


wave  coordinates  are  not  available  in  this  report.  Cleary  [17,  1965] 
measured  pressures  along  several  blunt  cones  at  M  =  5.25,  7.4,  and 
10.6.  In  addition,  pitot  pressures  were  .neasured  in  the  shock  layer 
for  a  15°  half  angle  blunt  cone  to  study  entropy  layer  thickness. 
Hasting,  Parsh  and  Redman  [37,  1957]  measured  body  surface 
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pressures  of  flat  faced  cones  and  provided  schlieren  and  shadow  photo¬ 
graphs  of  bow  waves. 

Experimental  body  pressures  were  compiled  by  Clark  [16,  1966] 
for  Mach  niunbers  from  1.  9  to  22  for  spherical  forebodies.  Empirically 
correlated  afterbody  surface  pressures  were  given  by  Eaves  [29,  1968] 
for  Mach  numbers  from  5  to  10.2.  Forebody  shapes  consisted  of  hemi¬ 
spheres,  flat  faces,  and  round  shouldered  flat  faces.  Bow  wave  coor¬ 
dinates  are  not  available  in  those  reports. 

Complete  experimental  boundary  layer  data  are  not  available  for 
axisymmetrical  blunt  bodies.  Boundary  layer  velocities  were  measured 
by  Wells  and  Blumer  [75,  1968]  for  a  hemisphere  at  M  =  2,  Reynolds 
number  per  inch  from  .  05x10^  to  .  5x10^  and  central  body  angles  of  30, 
50,  70  and  90  degrees.  Total  pressures  were  obtained  with  a  pitot  tube 
that  /as  normal  to  the  body  surface.  The  reported  uncertainty  of  the 
parameter  yD  was  ±  3  percent.  Here,  y  is  the  distance 

normal  to  the  body  surface,  D  is  the  maximum  diameter  of  the  body  and 
is  a  Reynolds  number  based  on  free  stream  static  conditions  and 
maximum  body  diameter.  Body  surface  cooling  was  provided  so  that 
surface  temperatures  approximated  the  wind  tunnel  stagnation  tempera¬ 
ture  to  approximate  an  adiabatic  body  surface  condition. 

2. 7  Summary 

For  inviscid  flow,  Aungier  [1,  1970]  was  successful  in  obtaining 
time  dependent  solutions  in  the  afterbody  region  of  a  blunt  body.  His 
results  agree  closely  with  measurements  and  those  of  the  method  of 
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characteristics.  However,  inviscid  solutions  of  blunt  body  systems 
cannot  provide  boundary  layer  and  heat  transft^r  information. 

Solutions  of  viscous  flow  about  blunt  bodies  have  been  limited  to 
two  dimensional  circular  cylinders.  In  the  case  of  Moretti  and  Salas 
[55,  1970],  solutions  were  restricted  to  the  forebody  region  with 
central  angles  less  than  70°  and  with  constant  viscosity.  Their  treat¬ 
ment  of  bow  waves  as  discontinuities  rcducod  the  computing  time.  The 
solutions  of  Scala  and  Gordon  [62,  1968]  required  excessive  computing 
times  and  were  not  verified  by  other  methods.  No  viscous  flow  solu¬ 
tions  were  reported  in  the  literature  for  the  afterbody  region  of  axi- 
symmetrical  bodies. 

Experimental  data  are  sparse  and  are  not  adequate  for  precise 
engineering.  The  most  complete  presentation  of  body  pressures  and 
bow  wave  coordinates  was  given  by  Baer  [4,  1961]  for  a  hemisphere 
cylinder.  For  blunt  bodies,  the  only  boundary  layer  experimental  data 
found  in  the  literature  were  those  of  Wells  and  Blumer  [75,  1968]. 


3,  THEORY 


3.  1  Introduction 

3. 1.  1  System.  For  purposes  of  this  investigation,  a  system  is  a 
specified  section  of  matter  that  is  being  considered  in  a  particular 
problem.  The  specifications  of  a  system  include  initial  values  and 
boundary  conditions. 

Specifications  for  supersonic  flow  about  a  blunt  body  system 
are  given  in  figure  1. 


Figure  1.  Specifications  of  the  Blunt  Body  System 
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In  developing  time  dependent  numerical  techniques,  Independent 

variables  -were  x,  r,  t,  and  dependent  variables  were  the  radius  of 

the  bow  wave  r^(x,t)  and  field  variables  fj(x,  r,t).  Boundary  conditions 

were  those  of  the  free  stream  and  body  surface.  Initial  values  r  (x,  0) 

w 

and  f.(x,  r,  0)  were  estimated  to  approximate  steady  values. 

3. 1.2  Scope.  The  problem  was  to  calculate  r^{x,  t)  and  f.(x,  r,t)  as 

t  -*  ®.  The  viscous  compressible  flow  equations  were  the  basis  of  cal¬ 
culations.  The  purpose  of  this  chapter  is  to  present  the  theory  thcii 
was  used  in  developing  numerical  techniques,  including  the  basic 
differential  equations,  coordinate  transformations  and  the  method  of 
characteristics. 

3.2  Basic  Differential  Equations 

From  [[64],  the  vector  form  of  the  basic,  viscous,  compressible 
flow  equations  are; 

If  +  V  •  (py)  =  0  (2) 

pV  =  B  +  V  •  S  (3) 

pTs  =  -v*Q  +  *:VV  (4) 

where:  p  is  density 
t  is  time 

V  is  velocity 

V  =  V  +  V  •  7V 

B  is  specific  body  force  per  unit  volume  which  is  negligible 
in  systems  of  this  investigation. 
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E  =  ”pl  +  #ls  stress  dyad 
p  =  pRT  is  pressure 
R  is  a  gas  constant 
T  is  temperature 
I  is  idemfactor 

Ss! 

♦  =  C(t?-  §M)V*  V]I  +2/4His  the  viscous  stress  dyad 

~  sa  fti 

17  is  bulk  viscosity  which  is  negligible  in  systems  of  this 
investigation,  [10] 

p  is  shearing  viscosity 

H  =  ^(7V  +  V7)isa  rate  of  strain  dyad 

V  7  is  the  transpose  of  TV 

s  is  entropy 

Q  =  -  k  7T  is  heat  flux 

k  is  thermal  conductivity. 

Conversion  from  vector  form  into  orthogonal  curvilinear  coor¬ 
dinates  was  accomplished  using  the  following  equations: 
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A 

h. 

i 


(5) 


TV 


hT  (i) 


(6) 


(7) 
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^  -r-^  -  6  yh  (8) 

dXp  dx^  a/3  a 

where:  repetition  of  the  English  indices  i  or  j  without  parenthesis 

implies  a  summation  of  that  term  with  indices  1,  2,  3 

h.  is  a  scale  factor  of  the  coordinates 
1 

y.  is  an  orthogonal  curvilinear  coordinate 

6  -  is  the  kronecker  delta. 

a/3 

3.  3  Coordinate  Transfornnations 

3.  3.  1  Introduction.  To  simplify  the  solution  of  the  blunt  body  system, 
body  related  coordinates  X  and  o)  were  used.  In  terms  of  them,  the 
entire  surface  of  any  body  is  specified  by  =  0.  In  the  numerical 
solution  of  equations  2,  3  and  4,  a  concentration  of  nodes  in  the  bound¬ 
ary  layer  and  near  the  stagnation  point  i^  desirable  for  accurate  resolu¬ 
tion  of  those  regions.  This  was  accomplished  with  nonlinear  coordinate 
transformations.  In  addition,  treatment  of  the  bow  wave  as  a  disconti¬ 
nuity  was  simplified  by  representing  the  bow  wave  as  a  grid  line. 

3.3.2  Body  Related  Coordinates.  The  body  related  coordinates  X 
and  <0  are  shown  in  figure  2,  where  X  is  the  distance  along  the  body 
from  its  nose  and  co  is  the  perpendicular  distance  from  the  body  surface 


to  a  point  in  the  system. 
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From  figure  2, 


sin  9  dX  +  (1}  cos  9 


cos  9  dX  -  (I)  sin  9 


(11) 


In  an  axisymmetrical  system,  velocity 

y  =  uX  +  v<3  (12) 

3.3.3  Nonlinear  Coordinate  Transformations.  In  the  numerical 
analysis  of  a  blunt  body  system,  greater  resolution  is  required  near 
the  stagnation  point  and  body  surface. 

Skoglund,  Cole  and  Staiano  [65]  and  Staiano  [67]  demonstrated 
that  nonlinear  coordinate  transformations  of  the  differential  equations 
yield  satisfactory  results  for  the  interaction  of  an  oblique  shock  wave 
and  laminar  boundary  layer.  The  logarithmic  transformations  of 
Skoglund,  Cole  and  Staiano  were  tried  for  the  blunt  body  system. 

They  were  later  modified,  so  that  the  bow  wave  was  represented  by  a 
grid  line.  This  simplified  coupling  of  the  bow  wave  to  the  digitized 
field. 

The  coordinate  transformations  that  were  used  in  the  solution 
phase  of  the  subject  investigation  are; 

C  =  Jtn(aX  +  1)  (13) 
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U  =  ln\$ - 7T~rr  +  0 

rw^(x,t)  j 


T  =  t 


where:  co  (X,t)  is  the  bow  wave  coordinate 
w 


a  and  ft  are  constants 


^i.t  ^  ^^i,A 

Xo) 

£,  =  (f.  ) 

i.  r  i,  r  - 
Zv 


(14) 

(15) 


For  equal  increments  of  AC  and  At/ ,  these  transforr  entrate 

nodes  near  the  stagnation  point  and  body  surface.  The  i.  ave  is 
represented  by  the  grid  line  v  =  ~  constant. 

3. 4  Transformation  of  Differential  Equations 

3.4. 1  Introduction,  The  basic  flow  equations  with  diffusion  were  not 
available  in  the  literature  in  terms  of  body  related  coordinates.  The 
simplest  method  of  expressing  equations  2,  3  and  4  in  body  related 
coordinates  is  to  use  the  scale  factors  of  equation  10.  A  more  tedious 
method  is  to  use  the  chain  rule  for  a  total  derivative.  As  a  check, 
both  methods  were  used.  In  these  derivations,  the  inviscid  and  diffu¬ 
sion  parts  were  treated  separately.  After  expressing  equations  2,  3 
and  4  in  terms  of  body  related  coordinates  \,  u)  and  t,  they  v/ere 
transformed  into  C>  V  and  r  using  equations  13,  14  and  15. 

3.4.2  Inviscid  Equations.  For  inviscid  flow,  ♦  =  0  and  Q  =  0.  In 

Si!  ^ 

transforming  equations  2  through  4  to  a  dimensionless  form,  reference 
variables  are  the  maximxim  body  radius  r^,  free  stream  acoustic 
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e  =  “^(m  sin  6  +  n  cos  0). 

0<P  pr 

The  symbol  e  is  used  because  this  term  appears  In  the  xate  of 
<P<P 

strain  tensor. 

In  cartesian  coordinates,  the  conservative  form  of  the  flow  equa¬ 
tions  Is  f,  ^  =  g,  +  h,  .  Lax  and  Wendroff  [45]  demonstrated  that 
l,t  I,  X  I,  y 

the  difference  form  of  that  equation  satisfies  the  Rankine-Hugoniot  re¬ 
lationships  for  a  finite  element  and  that  truncation  errors  of  these  dif¬ 
ference  equations  are  dissipative.  Equation  17  is  in  conservative  form 
except  for  the  term. 

From  equations  13  through  15, 


f 

I  V 


(18) 


V,c 


k-f 
3  ,v 


(19) 


f 

.  60 


4  ,  V 


(20) 
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(21) 

(22) 

^23) 


(24) 


From  equations  17  through  24 


f,  =  k,f,  +T^G,  . 
I,  r  1  I,  V  h  I,  C 


—G,  +  k  H,  +  B, 

h  1, 1/  4  I,  u  i 


3.4.  3  Viscous  Equations.  Viscous  flow  equations  2  through  4  have 
the  form 


Ic 

=  ^1^  +  irC.  .  -  -rrC.  +  k.H,  +  B,  +  v>, 

i,  T  1  i.  1/  h  i,  C  hi,!/  4  i,v  i  ’^i 


where  =  0 

+  ukp-  =  7  •  ♦ 

9  .  =  -7  •  Q  +  ♦  ;  7V 
4  ~  ~ 

The  derivation  of  these  terms  in  body  related  coordinates  was  lengthy. 
A  few  examples  of  the  procedures  are  presented  here.  From  equa* 
tion  6, 
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From  equations  7  and  8, 
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7  .  (TS)4[T(e,_^\.e^^<3)]_, 


+  rT(e  A  +  e  ,(S)]  +-{Te  0) 

wX  WO)  ,  (0  r'  ifitfr  ,<f> 


1  ^ 

*  lh'^'XX>,X^  — h  ^<^'xco>.7. 


+  — (e.  .  sin  0  +  •.  cos  0  -  e  sin  0)1  X 
r  '  XX  Xw  (p«p  M 


r<T6.  )  .  +  (T*  )  +-~(e 

l^h'  Xo)  ,  X  coco  ,0)  h 


coco  XX 


+  --(e.  sin  0  +  e  cos  0  -  e  cos  0)  ct)  (32) 
r  Xco  COCO  tp<p  ' 


By  similar  methods,  the  diffusion  terms  of  equation  26  are: 


<Pj  =  0 


M  C  , ,  1  .  _  _ 

*^2  ^  Re^  lh’-3^^®XX^X  ”  3^^•coco^X  ‘  3^^®(p<p\x 


+  4T*.  x]  +  2(Te-  )  +  — (3ln0e.  .  +  cos  0  e. 

Xco  Xco,  CO  r'  XX  Xco 


sin  0  e  )  I 
«}<p  I 


M  C 


^^3  “  Re  i'3^'^®Xx\co  3^'^*«co^co  "  3^’^®<p<p^co 


+  |(Te-  )  -  +^T(e  -  J  +  —  (sin  0  e, 

h'  Xoj  ,X  h  '  coco  XX  r  '  Xco 


+  cos  0  e  -  cos  0  e  )} 
coco  (pp  ! 


Tjafij-,  ';>■!?.  , -.  ■  ^  ,-v -;^»r 


M  C  , 

TRe  !pr-^^^^,X^X 
o  n 


+  T  cos  0)]  +  y(y  -  l)TC-f  (e^.  +  ef  +  e^ 

,  (0  3  A  A.  (0(0  <f><p 


-e»»e  -e-.e  -e  e  )+4e.  ] 

XX  (0(0  XX  <p<p  (0(0  <p(f>  X(0 


where:  M  =  free  stream  Mach  number 
o 

p  V  r 

Re  =  ■  —  is  free  stream  Reynolds  number 

o  u 
^o 

Pr  =  —r^  is  Prandtl  number 

K 

„  =.  C^T 

=  1  in  this  investigation. 

From  Sutherland’s  equation  [63], 


T  +  X  /T, 

C  -  8  su  I  b 

/  ~  T,  +  X  “IT 
b  su  »  s 

X  =  Sutherland  constant  =  198.  6®  R  for  air 
su 

T  =  stagnation  temperature 
s 

=  body  temperature 

In  transforming  equations  33  through  36  into  coordinates  Z 
V  of  equations  13  and  14,  the  second  derivatives  are: 


)  -  2h,k-f,  -  -  k_f, 

i,XX  3'  i,pi;  i,  p  2  3  i,  Cv  5  i,  i/ 


j.  r.  XV.3. 
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^l,  COW  ”  vv  "  ^i.  tf 


(38) 


f,  %  =  k,k A  .  •  -SfrC#  )  +  lc,k.(£.  -  £,  )  (39) 

l»COX  2  4\i.Cl/  ««.  £.w  3  4'l,vi;  i,v'  '  ' 


where:  ^2^  ^3*  ^4  given  in  cquettons  21  through  24. 


(40) 


w 


Complete  dU£erenee  equatlone  were  not  derived.  It  wet  simpler  to  uee 
equetione  18  through  40  directly  in  computer  codes. 

3,  4.4  Equations  Along  the  Stagnation  Streamline.  Equation  26  is 
Indetermlnant  along  the  stagnation  streamline  where  X  «  0  because  B| 
and  <p^  approach  •  as  r  approaches  0.  For  symmetry  about  the  stag¬ 
nation  streamline, 

'i'"-  'i.X  ■*'  'i.XX  W*®'  “'•»•«««> 


For  i  s  2, 


u  =  0,  to.  .0. 


(42) 


L'Hospital's  rule  is  that  i£  £|(X)  and  ^(X)  •*  0  as  X  **  0,  then 


limitl 

X^ 


In  the  £ollowing  example  o£  the  application  o£  I  'Hospital's  rule,  at 

j 

^  ^  «  “(u  sin  •  +  V  cos  f).  At  X  «  0,  r  «  C,  u  «  0,  0  ■  »/2  and 

*pv  * 
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Other  terms  of  equation  26  do  not  require  special  treatment. 

3.  5  Method  of  Characteristics 

Since  the  system  is  time  dependent,  the  bow  wave  moves.  In 
determining  bow  wave  speed,  one  equation  is  needed  in  addition  to  the 
Rankine-Hugoniot  relationships.  The  additional  equation  was  obtained 
by  the  method  of  characteristics.  In  the  following,  the  technique  of 
Moretti  and  Abbett  [52]  is  adapted  to  body  related  coordinates. 

Let  (i)  (X,t  +  At}  be  the  point  where  the  U)  coordinate  intersects 
the  bow  wave  at  time  t  -f  At.  In  this  section,  the  orthogonal  coordinate 
frame  of  figure  3  with  origin  at  <i)^(X,t  +  At)  is  used. 
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As  shown  in  figure  4,  the  boundaries  of  the  blunt  body  system  are 

A  B  C  D  E  F  A«  As  indicated,  the  coordinate  co  of  the  bow  wave 

w 

varies  with  respect  to  X  and  t.  Equation  14  implies  that  the  bow  wave 
coordinate  is  a  constant  and  is  not  a  function  of  ^  or  r.  In  the 
digitized  field  of  figure  5,  =  constant  is  one  of  the  boundaries. 

Transformations  of  cylindrical  coordinates  to  (  and  i>  were  specified 
in  equations  11,  13  and  14.  The  determination  of  a  and  of  equations 
13  and  14  is  described  in  section  4.2.2  because  of  their  dependence  on 
initial  values. 


4. 2  Initial  Values 

4.2. 1  Initial  Bow  Wave  Coordinates.  Initial  values  of  bow  wave 
coordinates  were  calculated  using  the  method  of  Love  [46],  The 
method  involves  a  combination  of  functional  and  empirical  analysis 
which  assumes  that  the  bow  wave  is  hyperbolic,  so  that 


r  (x,0)  = 
w 


nl/2 


M  -  1 
o 


(64) 


where;  (  =  distance  from  the  most  forward  point  on  the  bow  wave  to 

& 

an  intercept  of  its  asymptote  with  the  x-axls 
=  distance  from  stagnation  point  to  an  intercept  of  the  bow 


wave  asymptote  with  the  x^axis 


M  =  free  stream  Mach  number, 
o 
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A  characteristic  variable  P(ri,t)  is  introduced  so  that 
.  df 


,  t  dS 


f  - 

,V  “ 


Substituting  equations  52  and  53  into  50  and  51, 

dv 


B  +  V  B  +  oV  5  ^  =  c 

P,tdfi  ^n^.ri  d/J  d/5  ‘'l 


dv  , 

n  — tt  +  v  B  fl  ^  =  C 


where;  <=i  =  ^ 


c-  =  ‘•v„v  ^ 

2  €  t/.| 


The  determinant  form  is 


py^  _ 

>v 

d^ 

py^fV 

ft  t)  f  1) 

dv 

_J2. 

d^ 

"2 

(52) 


(53) 


(54) 


(55) 


(56) 


In  accordance  with  the  usual  method  of  characteristics  [53],  the  left 
determinant  is  set  equal  to  zero.  The  result  Is 

/3  +  (V  ±  a)/5  =  0  (57) 

ft  ff  ,rf 

One  solution  of  equation  57  is  the  characteristic 

(58) 


^  =  n  -  (v^  -  a)t 


I 

4 
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For  /3  constant, 


dn.  _ 

*  V  ■''»■* 


From  equations  50,  51  and  59  a  compatibility  equation  is 


,  dv 


(60) 


where  a  =  -pyv  -  v  p  +  ypav.v  . 
»)  §*5  %  i)»  % 


In  chapter  4,  equations  117  through  127  are  the  basis  of  a  wave  fitting 
technique  which  couples  the  bow  wave  to  the  digitized  field. 

In  this  investigation,  the  method  of  characteristics  was  also 
used  to  derive  an  inviscid  boundary  equation  at  the  body  surface  in 
terms  of  body  related  coordinates.  With  (o  normal  to  the  body  sur¬ 
face,  (t)  replaces  ij  and  v  replaces  v^  in  equations  59  and  60  so  that 


do). 


(61) 


dp.  dv, 

It’ff  ■  = »« 


(62) 


2 

where  =  rp<v  ^  '  h^.X  '  ’’'"'XX  ^  *w> 


The  value  of  in  equation  62  it.  indeterminant  along  the  stagnation 
streamline.  Using  the  technique  of  section  3.4.4,  at  X  =  0 


^(O^-^^'P^XX 


(63) 
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4.  NUMERICAL  TECHNIQUES 

4.  1  Introduction 

In  this  investigation,  the  transformed,  viscous,  compressible 
flow  equations  of  section  3.4  were  solved  numerically  by  converting 
them  into  a  set  of  difference  equations.  The  solution  was  started  with 
initial  values  at  all  nodes  and  Involved  boundary  values,  equations  for 
interior  nodes  and  equations  along  the  stagnation  streamline.  Those 
equation's  are  presented  in  this  section.  The  sequence  and  relationship 
of  computer  operations  are  presented  in  section  5. 

The  following  terms  are  used  in  the  description  of  numerical 
techniques.  A  system  is  the  specified  section  of  matter  that  is  being 
considered  in  a  particular  problem.  The  specifications  of  a  system 
include  initial  values  and  boundary  conditions.  A  field  is  the  Interior 
part  of  a  system  at  a  specified  time.  It  does  not  include  boundaries. 
The  dimensionless,  independent  variables  are  u  and  r  as  defined 
in  equations  13  through  15.  Nodes  are  defined  by  specified  values  of 
5  and  v.  The  Increments  AC  and  At/  between  nodes  are  constant. 
Intervals  are  changes  in  time  Ar  which  are  approximately  constant. 

A  cycle  is  a  set  of  computations  for  all  nodes  at  a  single  time.  An 
Iteration  is  one  step  of  a  successive  approximation  at  a  single  node  and 
time. 
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As  shown  in  figure  4,  the  boundaries  of  the  blunt  body  system  are 

ABCDEFA,  As  indicated,  the  coordinate  to  of  the  bow  wave 

w 

varies  with  respect  to  \  and  t.  Equation  14  Implies  that  the  bow  wave 

coordinate  v  is  a  constant  and  is  not  a  function  of  f  or  r.  In  the 
w 

digitized  field  of  figure  5,  =  constant  is  one  of  the  boundaries. 

Transformations  of  cylindrical  coordinates  to  (  and  u  were  specified 
in  equations  11,  13  and  14.  The  determination  of  a  and  jS  of  equations 
13  and  14  is  described  in  section  4.2.2  because  of  their  dependence  on 
initial  values. 

4. 2  Initial  Values 

4. 2. 1  Initial  Bow  W  ve  Coordinates.  Initial  values  of  bow  wave 
coordinates  were  calculated  using  the  method  of  Love  [46].  The 
method  involves  a  combination  of  functional  and  empirical  analysis 
which  assumes  that  the  bow  wave  is  hyperbolic,  so  that 


where:  =  distance  from  the  most  forward  point  on  the  bow  wave  to 

an  intercept  of  its  asymptote  with  the  x-axis 

=  distance  from  stagnation  point  to  an  intercept  of  the  bow 

wave  asymptote  with  the  x-axis 

M  =  free  stream  Mach  number, 
o 
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Love  empirically  determined  the  parameters  |  and  ^  and  tabulated 

3l 

them  ?s  functions  of  body  shape  and  free  stream  Mach  number. 

4.2.2  Coordinates  of  Nodes.  Equal  increments  of  AC  Al/ were 
used  in  computations  to  establish  the  location  of  nodes.  The  relation¬ 
ships  between  C>  T  and  indices  I,  J,  K  are: 


C  =  (I  -  DAC 

(65) 

V  =  (J  -  l)Ay 

(66) 

’’k+1  ^ 

(67) 

Cylindrical  coordinates  x(I,  J)  and  r{I,  J)  of  nodes  were  determined 
from  equations  11,  13  and  14. 

In  starting  a  viscous  flow  calculation,  the  number  of  nodes  in  the 
boundary  layer  and  the  total  number  of  nodes  between  the  body  surface 
and  bow  wave  were  selected.  Staiano  [67]  found  that  at  least  six  nodes 
were  required  in  the  boundary  layer  for  acceptable  accuracy.  The 
minimum  boundary  layer  thickness  is  at  the  stagnation  point.  From 
’.''I  the  boundary  layer  displacement  thickness  at  the  stagnation 
•u  .nt  axisymmetrical  flow  about  a  blunt  body  is 


where:  A.(g.)  =  the  tabulated  function  of  [18] 
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par 

Re  =-£-2-2 

sub'cript  8  Is  for  free  stream  stagnation  condt.  ?ons 

subscript  b  Is  for  the  body  surface 

subscript  e  is  for  the  edge  of  the  boundary  layer. 

Using  experimental  data,  Boison  and  Curtiss  [9]  derived  an  empirical 
equation  for  the  stagnation  point  velocity  gradient  u  .  .  Combining  it 

a  I A 

with  equation  68,  the  stagnation  point  boundary  layer  displacement 
thickness  is 


where:  =  value  of  X  for  p  /p  =  .  95 

e  s 

p  =  stagnation  point  pressure 
s 

From  Schlichting  [63], 

Wg  ^  36*  (70) 


where  is  velocity  thickness  of  the  boundary  layer. 

Using  the  Newton-Raphson  iteration  technique,  the  parameter  ft  of 
equation  14  was  calculated  from 

/  “ 

ej£n(/S  +  1)  -  1 

'  w 


0 


(71) 


where:  =  number  of  increments  in  the  boundary  layer  at  the  stagna¬ 

tion  po'-vt 

e_  =  number  of  increments  for  the  stagnation  streamline 

d  =  standoff  distance  of  bow  wave, 
w 

In  a  similar  manner,  the  parameter  a  was  calculated  using 

e.£n(0{X  +  1)  -  e  ln(aX.  +  1)  =  0  (72) 

J  ^  i  3 

where:  X  .  =  value  of  X  at  the  junction  of  the  forebody  and  afterbody 
J 

X  *  =  value  of  X  at  the  downstream  boundary 
I 

e^  =  number  of  increments  for  the  forebody 

e^  =  number  of  increments  for  the  entire  body  surface. 

4.2.3  Boundary  Layer.  A  simple  method  was  not  available  in  the 
literature  for  calculating  initial  values  in  the  boundary  layer  of  an 
axlsymmetrical  body.  From  an  extensive  review  of  the  liter^-'.ture, 
the  two  dimensional,  similarity  solution  of  Cohen  and  Reshotko  [19] 
was  selected  as  the  best  available  'T-asis  for  calculating  initial  values 
in  the  boundary  layer.  The  transformation  of  Mangier  [48]  was  used 
to  adapt  the  two  dimensional  solution  to  axlsymmetrical  flows.  The 
necessary  boundary  values  wer-.,  u  ,  p  ,  p  and  T  .  The  pressure 

6  G  G  G 

p  (X)  was  approximated  using  the  results  of  Clark  [16].  Assuming 

G 

isentropic  flow  along  the  outer  edge  of  the  boundary  layer,  u  (X), 

e 

p  (X)  and  T  (X)  were  easily  calculated  for  each  value  of  X(I).  The 

G  G 

results  of  Cohen  and  Reshotko  [19]  are  tabulated  in  terms  of  a 
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’^  =  '=«<S0s>'7  SPe-^^^ 

•^O 

f^T?^)  =  .9995 

M  =  Mach  number  at  o)  . 
e  e 

Equation  79  is  a  combination  of  transformations  of  Mangier  [48]  and 
Stewartson  [68],  In  the  computations,  derivatives  were  approximated 
by  centered  differences,  and  integrals  were  approximated  by  Simpson’s 
rule.  Values  at  nodes  were  determined  by  linear  interpolation  with 
respect  to  o). 

4.2,4  Initial  Field  Values.  The  initial  bow  wave  coordinates  were  cal¬ 
culated  using  equation  64.  The  Rankine-Hugoniot  relationships  were 
used  to  compute  initial  values  of  f^(X»  Ct)^>  0),  assuming  wave  speed 
equaled  zero.  This  provided  initial  values  for  fj(X ,  (0^»  0).  Utilizing 
initial  values  at  and  <o^,  by  linear  interpolation 

fj(X .  0),  0)  =  fj(X ,  (Og,  0)  +  C^[fj(X .  (0^,  0)  -  fj(X ,  cOg,  0)]  (80) 

where  C  =  (o)  -  CO  )/(co  -  to  ). 

to  ewe 

For  viscous  flows  tO  was  at  the  edge  of  the  boundary  layer.  For 

inviscid  flow  to  =0. 

c 

4. 3  Field  Equations 

4.  3.  1  Introduction.  Implicit  and  explicit  numerical  solutions  of 
equations  2  through  4  are  often  complicated  by  instabilities.  Many 
techniques  for  solving  them  are  available  in  the  literature  [1,  5,  8, 
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— . 


where:  is  an  arbitrary  stabilizing  term  that  is  described  in 

section  4.  5 

(C*  V,  t)  is  given  by  equation  26 

f,  (Z>  Vt  T)  is  derived  in  section  4.  3,  2 
l*  TT 

Depending  upon  the  particular  technique,  some  of  the  terms  of  equaticiik 
81  may  be  zero. 

4.3.2  Lax-Wendroff  Technique.  On  the  basis  of  previous  results  C54»' 
65,  66]  and  those  of  this  investigation,  an  extension  of  the  Lax- 
Wendroff  differencing  technique  [45]  was  used  to  solve  viscous  com¬ 
pressible  flow  around  a  system  with  an  afterbody.  The  main  feature 
of  the  technique  is  the  approximation  of  f^  in  equation  81.  Since 
V  =  u(X,  0),  t), 
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By  expanding 


0 

-1 

0 

0 

2  T 

u  +y(s-y) 

-2u 

0 

_T 

’y 

(84) 

uv 

-V 

-u 

0 

su 

-s 

0 

-u 

0 

0 

-1 

0 

uv 

-V 

-u 

0 

0 

-2v 

T 

"y 

(85) 

sv 

0 

-s 

-V 

By  expanding  equation  83, 


i,  tt  h  ij,X  ij,  <0  Ij  j,t 


a,. 


/tG-  %%  -“TG.  .h  .  +  H.  .  +  B.  a 

yh  j«  X  X  j^2  j»  X  >  X  j$  (oX  J’  Xy 


b.YrG.  .  --—G.  .h  +H.  +B.  \ 

ijlh  j.Xo)  j^2  j,X  ,co  J.WCO  j,  (Oi 


The  jacobians  a..,  b. .  and  c,.  were  differentiated  explicitly  to  avoid 

j  y 

storing  them  during  computations.  For  example,  from  equations  19 


and  84 


Si,  X  °  ’'V,  C  ■  V.  ^  “<V.  C  ■  V,  v’ 


An  example  of  the  differentiation  of  G,  H  and  B  is  that 


°3,X  ""  ^2^3,  C  "  NS,i/ 


4.3.3  Aungier  Technique.  The  differencing  technique  of  Aungier  [1] 

was  used  to  check  the  derivation  and  computer  coding  of  boundary  and 

field  equations.  The  technique  is  a  variation  of  one  proposed  by  Lax 

[44]  and  is  of  first  order  accuracy.  In  it,  f.  =0  and  0  in 

TT  I 

equation  82. 

4.3.4  Approximations  of  Derivatives.  Approximations  of  the  deriva¬ 
tives  that  evolve  from  equation  81  were  those  of  equations  18  through 
24,  37  through  40  and  the  following: 
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(88) 


i^j.{K.v,r)  =^Cfi(C+AC,v,T)-fj(C  -AC.y.T)] 

h,i;  " 

fi^j^(C»y.T)  =-^[VC  +  AC,i^,T)-2fj(?,i;,T) 

+  f.{C  -  AC.i/.t)]  (90) 

=  -^Cf,(C.  V  +  Av.t)  -  2£  (C,v,t) 

At/ 

+  f.(C.t/  -  Ai/.t)]  (91) 

fi_j^(c.x,T)  =:f^^aj(c  +  ac.i/  +  av,t)- 

fj(C  +  AC,  V  -  AI/,  t)  -  f.(C  -  AC .  +  Ai/,  t)  + 
fj(C-  AC.i/-  Ai/.t)]  (92) 


"w,  t  '  it  ‘  ^  ■  "w<^  '  <’’> 

"w.tt  =  iff"w,t<^-‘  +  "  -  "w,t<^-‘>’ 

ci0 

Because  X  =  is  discontinuous  at  the  junction  of  a  hemisphere- 

oA, 

cylinder,  the  approximations  of  derivatives  at  the  junction  were  not 
compatible  with  those  at  nearby  nodes.  To  avoid  this  problem,  an 
average  curvature 
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(95) 


xja  =  |cx(c  -  2AC)  +  x(C  -  AO  +  4X(C)  +  x{C  +  AC) 
+  X(C  +  2AC)] 

was  used  at  all  nodes. 

As  illustrated  by  figure  6,  there  is  no  discontinuity  in  x  . 

d> 


Coordinate  \ 

Figure  6.  Curvature  x  of  a  Hemisphere-Cylinder 
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(^i  xj.)»  =  (ira,.  .  +  b,.  +  c..)f. 

i,  tt'Xo)  'h  ij,X  ij.w  ij  j,t 

\  V  fr^,  -“tg.  xh 

p  h  I  X  c  li  yh  j  >  X  X  j^2  j » X  » I 

+  H.  »  +  B.  +b.Y^G.  . 

j.wX  j,Xy  ij^h  j.Xo) 

-^G.  .  h  +  H.  +  B.  \ 

J^2  J,X  ,  CO  J.COCO  J,  col 


^  \  t.  ~  “ZTZP  -nu.  -  3pu  .  e.  . 

*Xt  yhjXX  i\>co  »Xco  fXXX 


,3xv  ,  . 

'X-^'xx) 


(100) 


“<f''‘’xx>  " 


J.  4 

'T  *  'xx’ 


The  jacobians  and  b^^  of  equations  84  and  85  are  unchanged.  The 

derivatives  f,  and  f.  were  approximated  with  the  centered  differ - 

i,  V  1,  vv 

ences  of  equations  89  and  91.  The  derivatives  f j  \  »  fi  x  \  f.  » 

IiAA  iiAC0 

were  approximated  using  the  symmetry  conditions  of  equations  41  and 
42,  so  that  for  all  variables  except  m  and  u, 


4  =  ^4  %  =0 

l,X  I,  Xco 


(101) 


i.XX 


-2  Cfj(AX  ,v,r)  -  fj(0,  V,  t)3 


(102) 
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where  AX  =  [exp  (AC  “  1)3 /a. 


For  u, 


“.X 


(103) 


",xx  =  ” 


(104) 


X CO  ”  2AXAi/ ^  .  V  -  Av,  T)] 


(105) 


Derivatives  of  m  have  the  same  form. 


4.  5  Stabilizing  Terms 


4,5. 1  Introduction,  Linear  stability  theory  is  extensive  and  compli¬ 
cated  and  is  beyond  the  scope  of  this  investigation.  For  linearized 
one  dimensional  equations,  Richtmyer  and  Morton  [59]  derived  the 
following  criterion  for  numerical  stability  of  the  Lax-Wendroff  differ¬ 
encing  technique  when  0^  =  0: 


(|vj4a)^<, 


(106) 


where:  v^  =  the  x  component  of  velocity 
a  =  sonic  speed 

Equation  106  is  the  Courant-Fredrlcks-Lewy  criterion  for  stability. 
For  linearized  two  dimensional  equations  in  which  Ax  =  Ay,  Richtmyer 
[58]  derived  the  criterion  that 


/  2  .  2  ^  ^  1 
/v  +  V  +  a  I  ^  i  — 
X  y  /Ax  n 


(107) 
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where  v^  =  the  y  component  of  velocity. 

Numv.rical  experiments  by  Emery  [30],  Gary  [34]  and  others  have 
demonstrated  that  the  stability  criterion  of  equation  107  is  necessary 
for  many  systems.  It  was  used  with  the  equality  sign  in  this  investiga¬ 
tion  to  determine  the  upper  limit 

Al  - -  ^  -  (108) 

yi(^JT7.a) 

where  Ax  =  minimum  of  AX  and  AtO. 

Many  investigations  [12,  30,  31,  43,  66]  have  shown  that  linear 
stability  criteria,  such  as  equation  107,  are  necessary  but  not  suffic¬ 
ient.  In  previous  numerical  solutions  of  blunt  body  systems  In  which 
equations  106  or  107  were  satisfied,  numerical  instabilities  have  oc¬ 
curred  near  the  stagnation  point,  sonic  line  and  downstream  boundary. 
In  some  of  these  cases,  stabilizing  terms  were  added  to  the  basic  dif¬ 
ference  equations.  A  stabilizing  term  is  an  artificial,  mathematical 
term  that  is  added  to  a  difference  equation  to  improve  numerical  stabil¬ 
ity.  Satisfactory  stabilizing  terms  must  not  introduce  unacceptable 
errors  in  the  results.  In  this  investigation,  stabilizing  terms  of 
Lapidus  [43]  and  Aungier  [1]  were  modified  and  used  in  num>.  '  ;■  I 
experiments  because  they  were  satisfactory  in  previous  solutions  of 
inviscid  flow  about  blunt  bodies. 

4.5.2  Lapidus  Stabilizing  Term.  Skoglund  and  Gay  [66]  demonstrated 
that  the  addition  of  the  Lapidus  stabilizing  term  had  a  negligible  effect 


in  the  boundary  layer  of  a  flat  plate.  For  two  dimensional  inviscid 
flow,  Lapidus  [43]  analytically  demonstrated  that  his  stabilizing  term 
was  conservative  and  was  of  third  order  accuracy.  By  expressing  the 
Lapidus  stabilizing  term  of  [43]  in  C  and  u  coordinates,  one  of  the 
stabilizing  terms  that  was  used  in  equation  81  is 

+  Ly  (C .  f .  r)[k/._  T)[V^_  I 


where:  f  -(C.  V,  t)  =  77- [f,(C  +  AC.  V.  t)  -  f.(C.  l'* 

J.  V  flit  j  j 


j(C.i/.T)  =-3;f  -  f.(C  -  AC.u.t)] 


Lij(C  +  V,  r)  =  Cj  |u(C  +  AC.  t)  -  u(C,  i/,  t)|6 


Lij(C.  V  +  “.  t)  =  C^IvCC,  1/  +  Av,  t)  "  v(C,  u, 

Cj  and  are  arbitrary  constants.  Acceptable  values  were 
determined  by  numerical  experiment 

is  the  kronecker  delta  =  lifi=j;  =  0if  i?fj 
4.  5.  3  Aungier  Stabllljiiing  Term.  The  derivation  of  the  Aungier  sta¬ 
bilizing  term  [1]  begins  with  the  following  linearized  flow  equation: 


(110) 


where  Ot2~  constant. 

By  substituting  a  solution  of  the  form 

f  =  f  +  f .  exp(bt  +  cx)  (111) 

o  0 

into  a  difference  form  of  equation  110  and  requiring  |exp(bAt)|  ^  1, 
the  limiting 


(112) 


v^here;  f  =  constant 
o 

f.  =  constant  «  f 
0  o 

b  and  c  are  constants, 

Aungier  reasoned  that  the  Courant-Fredricke-Lewy  stability  criterion 
should  be  satisfied.  From  equations  106  and  110,  =  u  +  a.  This 

yielded 


Ax  u  +  a 


(113) 


For  two  dimensional  systems,  Aungier  assumed  a  stabilizing  term  of 
the  form 
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By  simple  replacement  for  invlscid  flow  In  the  subject  investigation,  a 
stabilizing  term  that  was  used  in  equation  81  was: 


^4  = 


C3(u 


+  ar 


ML 

2 


i.XX 


+  C^(v  +  a) 


^f 

2  i,  COCO 


(115) 


where  and  are  arbitrary  constants. 

This  equation  was  satisfactory  for  inviscid  flow  but  yielded  inaccurate 
results  for  viscous  flow.  An  accurate  solution  for  viscov.s  flow  was 
obtained  using 


CgCu 


ba)^  ^f.  .  .  +C,(v  +  ba)^ 
2  i,XX  o 


M. 

2 


i,  COb> 


(116) 


where  b  =  exp  -  J)At/]  for  J  ^ 

=  1  for  J  >  J 

e 

is  an  arbitrary  constant 
J  =  xndex  of  .loies  in  the  co  direction 
J  is  a  node  near  the  edge  of  the  boundary  layer. 

4.6  Bow  Wave  Techniqu^i 

4.  6. 1  Bow  Wave  Equations.  At  the  upstream  boundary  A  B  of  figure  4, 
the  free  str>'  r.m  values  are  constant.  Since  the  governing  equations  are 
time  depend ^.it,  the  bow  wave  may  move.  In  figure  3  of  section  3.5, 
the  speed  of  the  bow  wave  in  the  %  direction  is  designate  d  as  w.  The 
Rankine>Hugoniot  relationships  for  a  moving  wave  are 


y  +  1 
2 


w)^  +  1 


-  w) 


+  w 


(117) 
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^C2  "  "^€1 


PZ  = 


(118) 

(119) 


Pz  = 


where; 


=  1 

+  -^[(v 
y  +  1  7) 

V,  f  U)  . 

-  1  I—wJl 

r)l 

r  L  h 

V  r 

=  ~|cos  0  - 

r  L 

cos  6  +  sin  6 
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(120) 


Vj  =  freestream  velocity  in  x  direction 


r  =  [(^)  ^ 


.1/2 


subscripts  1  and  2  refer  to  values  at  upstream  and  down¬ 
stream  sides  of  the  bov.=  wave. 

Unit  vectors  are; 


fi  =  -<s)  /r 


(121) 


(122) 


In  the  above  set  of  equ"  'ons  to  .  was  approximat-  1  by 

A 
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CO  »  (\ , t  +  At)  =  TCto  »(X,t  +  At)  +  2co  ^  (X , t  +  At)] 

w,  A  j  w,  A  'f/i  A 


•^CcoJC  +  AC.T  + At) 


(123) 


+  to^(C.  T  +  At)  -  ZuiJZ  -  AC ,  r  +  At)] 


where:  co^  ^^(X.t  +  At)  =  ■^Cw^(C  +  AC.T  +  At)  -  w^(C.T  +  At)] 


“w.X^^'"  ^  "  “w^^  "  AC,t  +  At)] 

Aungier  [l]  found  that  the  technique  of  equation  123  inhibited  oscilla¬ 
tions  of  (0  with  respect  to  X* 
w 

Equation  59  was  used  to  obtain  the  characteristic  slope  of  figure  7, 


t  +  At 


Characteristic  /3  •  Const. 


dr»  . 

~  s  V  -  a 
dt  rf 


Coordinate  rj 

Figure  7.  Characteristic  at  the  Bow  Wave 


The  intersection  of  the  characteristic  with  the  rj  axis  at  time  t  is 
symbolized  by  rj*  and  was  calculated  by  a  method  that  is  outlined  in 
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the  next  section.  Along  the  characteristic,  an  approximation  of  equa¬ 
tion  60  is 

v^2  "  ■*■  ^ 

where:  ypa  =  .5[yp^&^+yp*aL*'i 

primes  refer  to  values  at  and  time  t 

subscript  2  refers  to  values  on  the  downstream  side  of  the 
bow  wave  at  time  t  +  At. 

In  evaluating  a* ,  known  values  at  nodes  were  interpolated  to  determine 
value  -c  tj”*.  Using  the  chain  rule  yielded. 


afi 

at 

=  ip" 

3  ■’’TP 

“4 

(125) 

where: 

= 

cos 

(9  -  a') 

"w,X 

h 

sin  (9  ’ 

■9')j 

/  r 

r 

- 

“4 

= 

sin 

(9-0') 

h 

cos  (9  ■ 

•  9')l 

/  r. 

At  TJ^ 

< 

ti 

u  'a.  - 
4 

V 

(126) 

u^a^  + 

(127) 

The  values  of  X  U)*  and  0^  at  rj*  were  determined  using  the  Newton- 
Raphson  iteratl.’n  technique. 

The  bow  wave  coordinates  at  time  t  +  t*"  were  calculated  using 
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f  *1-  rv 


co^{X,t  +  At)  =  w^(X,t)  -  rw(t)At 


(128) 


and  held  constant  during  a  cycle. 

4.6.2  Iteration  Techniques.  To  calculate  f.^CX.tO^it  +  At),  equations 
of  the  preceding  section  were  solved  by  an  iteration  technique.  As  a 
first  approximation  of  w(t  +  At),  w(l)  =  w(t)  where  the  n^mber  In 
parentheses  refers  to  the  approximation  number.  Having  the  approxi¬ 
mate  value  of  w(l),  equations  17  7,  119  and  120  were  solved  for  v  ,(1), 

tjZ 

p^Cl)  and  p2(l)*  With  available,  the  first  approximation  of  f)* 

was 

t?^(l)  =  -  ^^(DlAt  (129) 

where  ^2^1)  =  yr^TT)* 


Using  the  Newton-Raphson  iteration  technique  and  linear  interpolation, 

a^  was  determined  and  equation  126  was  solved  for  v* .  The  second 

V 

approximation  of  rf*  was 


n'(2)  =  \v^^(l)-a^(l)+v^-a'\At/2 


(130) 


The  above  process  was  repeated  until 

|t7'(n)  -  Tj^n  -  1)|  <  10  ^ 


(131) 


Once  an  accurate  location  of  n*  was  obtained,  o*  was  calculated  and 

V 

equation  124  was  solved  for  v  «(2). 

TJZ 

Ti 


65 


. - 


>10'^  (132) 

a  new  value  of  w(2)  was  estimated  from 

»^2)  +  C^v^j  -^v^j,(l)W2) 


The  value  of  w(2)  was  used  in  equations  117,  119  and  120  to  estimate 

V  ~(2),  p_(2)  and  pA2).  A  new  value  of  corresponding  to  v  -(2)  and 

1J2  2  2  *  ty2 

a^CZ)  was  determined  from  equations  129  through  131.  Equation  124 

was  again  solved  for  v  .(3).  This  process  was  continued  until 

rfi 


(134) 


Only  a  few  iterations  were  required  to  satisfy  equation  134. 
After  satisfying  equation  134,  the  downstream  X  and  co  components  of 
velocity  were  computed  from 


«(x,«„.t4/)t).[v^2^4vjJ/r  (13S) 

v(X.»„.l  +  /«)  =  (136) 

4.  7  Boundary  Equations 

4.  7. 1  At  the  Downstream  Boundary.  The  difference  equations  at  the 


downstream  boundary  X  ^  were  the  same  as  those  in  the  field  of  section 
4.3.  Values  at  X.,.  were  required  in  the  approximations  of  derivatives 
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3-t  In  finding  the  location  of  nodes  at  an  approximation  for 

shock  wave  coordinate  U3  (X..,)  was  needed.  From  figure  8,  the 
extrapolated  value  of  was  determined  by  finding  the  intersec¬ 

tion  of  a  line  drawn  through  the  points  u)  (X ,  , )  and  03  (X  ,)  with  a  line 

W  A  w  X 

extended  in  the  (O  direction  at  X,.  ,♦ 

X I  A 


Bow  Wave 


Figure  8.  Location  of  co  (X  .  , ) 

w  i+1 


From  the  law  of  sines, 

Ww-t^jj^.,)  =  b  sin  a  /  sin  /S  (137) 

By  linear  extrapolation, 

=  V^^,j.w.t)  +  C^Cfj(X^,w,t)  -  fj(X^_j,a),t)] 

(138) 
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whei*e  "  X^_j) 


4.  7.  2  At  the  Body  Surface  with  Viscous  Flow.  The  boundary  condi¬ 
tions  for  viscous  flow  at  the  body  surface  were  u(X ,  0,t)  =  v(Xi  0,  t)  =  0 
and  T(X,  0,t)  =  =  constant.  Using  s  =  jtn(p/p')  to  combine  equations 

for  i  =  1  and  4  in  equation  26, 


".t  “ 


(139) 


where 


y  r  2 

:  <p.  =  +T.T 


+  XT^T 
(t)U}  b  ,  0) 


60  1 
—  COS  0  I 


+  y{y  -  +u^  ] 

b  3  60(0  ,  (0 


At  the  stagnation  point 


M  C  (  - 

<p.  =  +T,T  +2xT,T  ] 

^4  Re  Pr  ,  60  b  ,  60(0  b  ,  60 


(MO) 


+  |y(y  - 

3  b  6060 


Using  a  Taylor  series,  the  pressure  at  the  body  surface  is 


p(X,0,t  +At)  =  p(X,0,t)  +  Atp  .(X.0,t) 

*  ^ 


(Ml) 


where:  p  ...  =  -yp  ,v  -  ypv 

tt  '  t  ,  60  ,  60t 


2  2 

V  .  =  -V  +  p  p  /yp  -  p  /yp 

,60t  ,60  ^,  60  ,  60  %60(O 


-  ^,,1;  r"‘-  . 


(p  =  <p  .  =  0  as  an  approximation 

^  ^  (a)  4 1 1 

(142) 

i«  (0  Aco  1  i 

hu)  =  W^Cexp(Au)  -  1]/^ 

p  ^^(X  I  0,  t)  was  approximated  by  equation  38  for  o)  =  Aco* 

The  density  at  the  body  surface 

p(X.0,t  +  At)  =  p(X,0.t  +  At)/Tj^  (143) 

4.  7.  3  At  the  Body  Surface  with  Inviscid  Flow.  Boundary  conditions  at 
the  body  surface  for  inviscid  flow  are  v(X,  0,  t)  =  0,  u(X>  0,  t)  ^  0  and 
T(X,0,t)  is  not  constant.  For  inviscid  flow, 


,t  =  -““.X 

(144) 

(145) 

Equation  6l  was  used  to  establish  the  characteristic  of  figure  9. 
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Figure  9.  Characteristic  at  Body  Surface 


The  intersection  (a*  on  the  (o  axis  at  time  t  was  found  using  the  iteration 
techniques  of  section  4.6.2.  Equation  62  was  approximated  by 

p(X,  0,t  +  At)  =  0"^^  +  p^  -  ypa  v'  (146) 


where:  ypa  =  .  5y  Cp{X  ,  0,  t)a(X  ,  0,t)  +  p'a '] 
primes  refer  to  values  at  . 

Equations  144  and  145  were  solved  in  the  same  manner  as  equation  141. 

The  derivative  f.  was  approximated  by  equation  142,  and  the  deriva- 
CO 

tive  f.  was  approximated  by  equation  38  for  co  =  Au).  Equations  19 
COCO 

and  37  with  k-  =  k_  =  0  were  used  to  calculate  f.  .  and  f.  ^  .  In  these 

3  5  i,X  i,XX 

equations,  the  derivatives  f.  ^  and  f.  ««  were  approximated  by  equa- 
and  90. 
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f 


At  the  stagnation  point,  u  =u  =s  =  s  =0.  For  invlscid 

f  W  f  CC  f  V  I  Vb 

flow  the  stagnation  point  temperature  was  precisely  calculated  from  an 
isentropic  equation  and  was  constant. 


5.  COMPUTATIONS  AND  RESULTS 


5.  1  Introduction 

In  this  investigation,  digital  computers  were  used  to  solve  the 
difference  equations  presented  in  section  4.  In  the  first  computer 
codes,  there  were  major  problems  with  nodal  distribution,  wave  fitting 
and  numerical  instabilities.  Refinement  of  difference  equations  and 
computer  codes  was  accomplished  experimentally  on  digital  computers. 
In  order  to  keep  the  length  of  the  report  within  a  reasonable  limit, 
only  the  more  significant  computations  and  results  are  presented.  The 
two  principal  computer  codes  are  INITIAL,  based  on  initial  value 
techniques  of  section  4.2,  and  MAIN,  based  on  the  techniques  of  sec¬ 
tions  4.  3  through  4.  7.  Results  are  presented  for  axisymmetrlcal, 
inviscid  and  viscous  flows  around  hemispheres  and  hemisphere- 
cylinders.  Because  of  stability  problems  and  constraints  on  computa¬ 
tional  time,  most  development  computations  were  restricted  to  hemis¬ 
pheres.  Table  i  lists  specifications  of  the  systems  that  were  solved 
in  this  investigation. 


Table  1.  Systems 


System 

Number 

Configuration* 

Reynolds 
Number,  R  e 

o 

Mach 

Number,  M 

o 

1 

K  -  C 

Inviscid 

4 

2 

H 

10^ 

4 

3 

H  -  C 

4  X  10^ 

4 

4 

H 

lo" 

4 

- 

H 

10^ 

4 

6 

H 

Invijcid 

2 

7 

H 

10^ 

2 

*  Under  "Configurati'^n,  "  H  =  hemisphere  a  me;  H  -  C  =  hemisphere- 
cylinder  combination. 

For  most  of  the  systems  of  table  J,  many  cor.xputer  runs  were 

required  before  results  ere  satisfactory.  Systems  1  and  4  were  used 

as  test  cases  to  develop  the  MAIN  code  and  stabilizing  terms.  Aungier 

[1]  solved  system  1  by  segmenting  the  field.  As  a  heck,  his  results 

were  compared  with  those  of  this  investigation.  Steady  or  time 

dependent  solutions  of  viscous  systems  with  an  afterbody,  such  as 

system  3,  were  not  previously  available  in  the  open  literature.  Inves- 

tigatiix  of  the  effect  of  Reynolds  number  at  M  =4  was  the  motivation 

o 

for  solving  systems  2  and  4.  Solutions  of  systems  2,  3,  4,  5  and  7 
demonrtrated  that  the  codes  were  satisfactory  for  viscous  systems. 
Solutions  for  systems  6  and  7  demonstrated  that  the  codes  are  satisfac¬ 
tory  for  lower  Mach  numbers  which  have  less  inherent  stability. 
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The  specifications  for  computer  runs  are  presented  in  table  2. 
Results  of  these  runs  are  presented  and  discusfied  in  succeeding  sec¬ 
tions.  In  table  2,  IJ  equals  I  just  ahead  of  th  *  junction  of  the  hemi¬ 
sphere  and  cylinder,  EL  is  at  the  downstream  boundary,  JE  is  at  the 
edge  of  the  boundary  layer,  JW  is  at  the  bow  wave  and  KL  is  when  com¬ 
putations  were  halted.  In  the  following,  most  figures  of  a  sub-subsec¬ 
tion,  such  as  5.  2. 1,  are  presented  at  the  end  of  that  section. 

5.2  Computer  Codes 

5.2. 1  INITIAL.  A  simplified  flow  diagram  of  the  INITIAL  code  is 
presented  in  figure  10.  In  it,  ILPl  =  IL+l  and  JEPl  =  JE+1.  Both 
magnetic  tape  and  punched  cards  were  used  to  input  data.  The  bound¬ 
ary  layer  tables  of  Cohen  and  Reshotko  [19],  bow  wave  parameters  of 
Love  [46]  and  body  surface  values  of  Clark  [16]  comprised  several 
thousand  entries  which  were  recorded  on  magnetic  tape.  Punched  card 
input  was  used  to  specify  body  shape,  reference  variables,  boundary 
values  and  nodal  parameters.  After  calculating  initial  values  of  p, 
m,  n  and  S  at  each  node,  an  output  tape  of  their  values  was  generated. 
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CF  »  Configuration:  H  »  hemisphere,  H-C  «  hemisphere-cylinder.  2.  Diff.  T  ■  Differencing  Technique. 


Execution  time  of  program  INITIAL  was  approximately  2  minutes  on  a 
CDC  3800  computer. 

Initial  bow  coordinates  were  computed  by  the  techniques  of  sec¬ 
tion  4. 2. 1.  The  computed  bow  wave  was  slightly  upstream  of  the  ex¬ 
perimental  one  of  Baer  [4]  as  displayed  in  figure  11. 


0  12  3  4 

Coordinate  x 

Figure  11.  Computed  and  Experimental  Bow  Wave 

Coordinates  at  M  =4 
o 

Since  the  MAIN  code  adjusted  the  location  of  the  bow  wave,  the 
difference  between  computed  and  experimental  bow  wave  coordinates 
was  not  important. 

Initial  values  of  pressure  in  the  boundary  layer  were  assumed  in¬ 
dependent  of  Reynolds  number  and  o).  As  expected,  computed  and 


T'"*'t'^T 'T' 


experimental  initial  body  pressures  of  figure  12  are  in  good  agreement, 


since  the  INITIAL  code  used  experimental  results  of  Clark  [16]  for  a 


hemisphere.  In  the  afterbody  region,  the  initial  pressure  Pjj(^') 


equal  to  the  initial  pressure  Pj^  at  the  junction  of  the  hemisphere  and 


cylinder.  As  shown  in  figure  13,  this  assumption  resulted  in  a  per¬ 


turbation  in  the  velocity  thickness  of  the  boundary  layer  near  the 


junction.  In  the  time  dependent  computations  of  the  MAIN  code,  these 


perturbations  in  co^  were  quickly  eliminated. 


ssau^iainx  i^Avi  Aa^punog 


Figure  13.  Initial  Values  of  Velocity  Thickness 
of  Boundary  Layer  for  M_  =  4  and  Re_  =  4000. 


5.2.2  MAIN,  A  simplified  flow  diagram  of  the  MAIN  code  is  pre¬ 
sented  in  figure  14.  A  set  of  indicators  was  used  to  con'-rol  input  from 
the  initial  value  tape,  specify  inviscld  or  viscous  flow,  specify  the  sta¬ 
bilizing  term,  set  the  output  and  run  intervals  and  specify  the  size  of 
the  field.  Integer  values  of  I,  J  and  K  were  used  to  identify  nodes  and 
time.  The  relationships  between  C*  t'*  T  and  I,  J,  K  are  given  by 
equations  66,  67  and  68.  KIN  and  KL  are  starting  and  ending  indices 
of  times.  If  KIN  =  1,  initial  values  were  obtained  from  INITIAL,  If 
KIN  >  1,  initial  values  were  obtained  from  an  output  tape  of  a  previous 
run.  Computations  were  halted  at  KL,  and  results  were  printed  and 
recorded  on  tape. 

In  operation  1  of  figure  14,  equation  109  was  used  to  compute 


At  -  at  all  nodes.  At  equaled  the  minimum  of  At ,  for  one  cycle  of  the 
A  loop.  In  operation  2,  bow  wave  coordinates  were  updated  according 
to  section  4,  6.  2.  Operation  3  computed  values  at  the  body  surface,  and 
operation  4  computed  values  along  the  stagnation  streamline.  In  oper¬ 
ation  5,  values  at  all  other  nodes  were  calculated  using  the  equations 
of  sections  4.  3  and  4.  5.  Extrapolation  of  values  at  the  downstream 
boundary  was  accomplished  in  operation  6. 
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2.  CALCULATE 
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SECTION  4.6 


4.  CALCULATE 
f.  (I.  J,  K) 


SECTION  4.4 


5,2.  3  Computational  Times.  Field  values  were  computed  at  a  rate 
of  56.  7  nodes /second  on  a  CDC  3800  computer.  Values  of  f^(I,  JW,  K) 
at  the  downstream  side  of  the  bow  wave  were  computed  at  a  rate  of 
17.8  nodes/second.  For  a  CDC  3800  computer, 

t^  =  (K)(IL)[c(JW  -  1)  +  .0563]  (147) 

where:  t  =  execution  time  in  seconds  ±  5% 
c 

c  =  .  0154  for  inviscid  flow  and  =  ,0177  for  viscous  flow 
K  i  100  cycles 

The  percentage  of  total  CDC  3800  computer  time  required  for  wave 


fitting  is  shown  in  figure  15. 


Number  of  Nodes  in  u)  Direction  JW 


Figure  15.  Percentage  of  Computer  Time  for  Wave  Fitting 
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In  this  investigation,  JW  =  11  for  inviscid  flows  and  16  a  JW  i  29  for 
viscous  flows.  Assuming  that  wave  fitting  could  be  accomplished  as 
quickly  as  field  computations,  from  10  to  48  percent  of  the  total  compu¬ 
tational  time  could  be  saved.  Numerical  experiments  were  made  in  an 

attempt  to  speed  up  the  wave  fitting  computations.  The  attempt  first 

-5  -3 

involved  increasing  10  to  10  in  expression  132.  The  second  at¬ 
tempt  was  to  use  equation  129  to  compute  rj*.  Neither  of  those  attempts 
was  successful  because  of  lack  of  convergence  near  the  bow  wave. 

5.  3  Stabilization  of  Computations 

5.3.1  Introduction.  The  first  computation  of  viscous,  compressible 
flow  around  a  hemisphere-cylinder  without  stabilizing  terms  displayed 
numerical  instabilities.  A  major  problem  of  this  investigation  was  to 
achieve  numerical  stability  without  destroying  accuracy.  Many  attempts 
were  made  to  avoid  the  use  of  a  stabilizing  term,  including  one  dimen¬ 
sional  wave  fitting  investigations,  stationary  coordinate  frames,  modi¬ 
fication  of  coordinates  of  nodes  and  variation  of  differencing  techniques. 
Results  of  unstabilized  computations,  the  development  of  satisfactory 
stabilizing  terms  and  accurate  solutions  of  several  systems  are  pre¬ 
sented  in  this  section. 

5.3.2  Unrtabilized  Computations.  Specifications  of  runs  1  through  4  for 
viscous  and  inviscid  flows  with  ^^  =  0  are  given  in  table  2.  The  location 
and  extent  of  instabilities  are  displayed  in  figures  16  through  19  at  the 
end  of  this  section  on  pages  87  to  90,  For  viscous  flow  in  run  1,  spatial 
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oscillations  in  the  u-velocities  at  I  =  3  of  figure  16  began  near  the  center 
of  the  field  and  spread  into  the  region  next  to  the  body.  As  displayed  in 
figure  17,  there  was  an  instability  in  u-vclocities  in  run  2  at  I  =  19. 
There  were  similar  instabilities  in  other  variables.  These  instabilities 
were  maximum  at  I  =  19.  Oscillations  were  small  near  the  body  surface 
and  were  maximum  near  the  center  of  the  field.  Results  of  runs  3  and  4 
for  inviscid  flow  are  shown  in  figures  I*"-  and  19.  Since  the  instabilities 
occurred  for  inviscid  flow,  they  were  attributed  to  the  inviscid  part  of 
the  field  equations.  Since  the  instabilities  were  similar  for  viscous  and 
inviscid  flows,  it  was  concluded  that  viscous  terms  had  little  effect  on 
the  instabilities. 
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u-veloclty 

Figure  16.  Instabilities  In  Forebody  Region  at  I  = 
Viscous  Flow  Around  a  Hemisphere  with  =  0. 
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Figure  17.  Instabilities  in  Afterbody  Region  at  I  =  19.  Run  2, 
Viscous  Flow  Around  a  Hemisphere- Cylinder  with  ^  =  0, 


u-veioclty 

Figure  It.  Instability  in  Forebody  Region  at  I  =  3. 

Run  3.  Inviscid  Flow  Around  a  Hemisphere  with  ^  =  0 
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5.  3.  3  Development  of  Stabilizing  Terms.  A  development  program  was 
initiated  to  determine  the  cause  of  the  instabilities  and  to  find  techniques 
for  stabilizing  solutions.  Inviscid  flow  solutions  converged  approximately 
four  times  faster  than  the  corresponding  viscous  ones.  Therefore, 
inviscid  systems  were  used  in  this  phase.  The  first  step  was  to  deter¬ 
mine  if  instabilities  were  caused  by  errors  in  equations  or  computer 
coding.  To  detect  errors  in  equations  of  section  4.  3.  2  and  their  coding, 
results  of  computations  using  them  were  checked  with  those  of  a  simpler 
two-step  technique  of  the  form 

+  =  fi(C.i/.r)+-^fj  ,j.(C.t/,T)  (148) 

fi(C.t/.T  + At)  =  fj(C,t/,T)  + Arfj  .j.(C.v.r+-^)  (149) 
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Using  these  equations  with  the  specifications  of  runs  3  and  4  of  table  2, 
the  instabilities  were  similar  to  those  of  figures  18  and  19,  except  they 
grew  more  rapidly.  It  was  concluded  that  the  equations  of  section  4.  3.  2 
were  correct  and  properly  coded. 

The  differencing  technique  of  section  4.  3.  3  with  a  stabilizing  term 
was  used  in  run  5  of  table  2  to  determine  if  wave  fitting  and  downstream 
boundary  equations  introduced  perturbations.  The  solution  was  stable 
at  K  =  800.  Computed  body  pressures  of  figure  20,  page  96,  are  in 
agreement  with  the  measurements  of  [4],  The  differences  of  computed 
and  measured  bow  wave  coordinates  in  figure  21  were  not  considered 
detrimental  to  other  parts  of  the  field  and  were  partially  attributed  to 
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the  averaging  of  body  curvature  x(X)  in  equation  96.  It  was  concluded 
that  the  bow  wave  and  downstream  boundary  equations  were  satisfac¬ 
tory.  The  results  of  runs  5  and  6  with  the  differencing  techniques  of 
sections  4.3.3  and  4.3,2  and  stabilizing  terms  were  within  1  percent 
of  each  other  except  near  the  junction  of  the  hemisphere  and  cylinder. 
The  results  of  run  6  were  smoother  and  in.  better  agreement  with  the 
measurements  of  [4]. 

Equation  110  was  used  to  stabilize  the  solutions  in  runs  7  and  8. 

In  run  7,  C,  =  C_  =  1,  and  in  run  8,  C.  =  C_  =  3.  In  figure  22,  body 
surface  pressures  of  those  runs  display  an  incorrect  hump  near  the 
junction.  Body  surface  pressures  of  run  6  in  figure  22  are  in  better 
agreement  with  the  measurements  of  [4]  than  are  those  of  runs  7  and  8. 
As  shown  in  figure  23  near  the  junction  of  the  hemisphere -cylinder, 
the  V- velocities  of  run  8  deviate  from  those  of  runs  6  and  7.  From 
figures  22  and  23  for  inviscid  flows,  the  stabilizing  term  of  equation 
115  yielded  more  accurate  results  than  those  of  equation  110, 

The  development  of  stabilizing  terms  for  viscous  flows  was  more 
difficult  than  for  inviscid  flows  because  conditions  in  the  boundary  layer 
are  sensitive  to  stabilizing  terms.  In  run  9,  with  the  differencing  tech¬ 
nique  of  section  4.  3.  3  and  th"  stabilizing  term  of  equation  115,  results 
were  unstable  near  the  stagnation  point  at  KL  =  68,  Run  10,  with  the 
differencing  technique  of  section  4.3.2  and  stabilizing  term  of  equation 
115,  was  stable  but  was  unsatisfactory  because  the  stabilizing  term 
was  much  larger  than  the  viscous  term  ip.  in  the  boundary  layer.  This 
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meant  the  boundary  layer  results  were  inaccurate.  A  few  runs  wore 
tried  with  a  new  stabilizing  term 


(150) 


where  =  C-  =  1,  2  or  3,  Although  0.  was  small  compared  to  <p.  in 

0  7  t  t 

the  boundary  layer,  severe  oscillations  occurred  outside  the  boundary 

layer  and  equation  150  was  abandoned. 

Bun  11,  for  viscous  flow  around  a  hemisphere,  utilized  the  stab¬ 
ilizing  term  of  equation  116  for  which  0^(X,  0,  t)  =  0  at  the  body  surface. 
As  shown  in  figures  24,  25  and  26,  in  the  boundary  layer,  the  stabilizing 
term  0.  is  small  compared  to  the  viscous  term  <p^.  Small  oscillations 
in  ^2  were  not  reproduced  in  figure  25,  The  steady  bow  wave  coordl- 
nat  5  and  body  pressures  agreed  with  the  measurements  of  [4].  The 
results  of  run  11  are  given  in  section  5.5.2.  In  runs  11,  and  17  through 

20,  the  stabilizing  term  of  equation  116  was  satisfactory  for  hemispheres 

3  5 

at  M  =2  and  4  and  Re  =10  to  10  ,  For  each  set  of  Mach  and 
o  o 

Reynolds  numbers,  numerical  experiments  were  used  to  determine  val¬ 
ues  of  Cg,  Cg,  C^,  and  nodal  parameters. 

For  viscous  flow  around  a  hemisphere-cylinder,  the  stability  prob¬ 
lem  was  more  difficult  than  for  a  hemisphere.  Figure  27  displays  pres¬ 
sures  versus  w/co^  in  the  afterbody  region  at  I  =  19  for  runs  2,  12  and 
13.  Run  2  was  unstabilized  and  was  discussed  in  section  5.  3.  2.  Runs 
12  and  13  used  the  differencing  technique  of  section  4.  3.2  and  stabilizing 
terms  of  equations  116  and  110  respectively.  Pressures  in  run  12  at 


KL  =  2000  are  unstable  and  are  approaching  those  of  run  1  without  stab¬ 
ilization.  In  run  13,  small  oscillations  are  present;  however,  these 
oscillations  decreased  as  K  increased.  Values  of  other  variables  were 
smooth  in  run  13.  As  illustrated  in  figure  28  for  run  13,  0,  was  larger 
than  in  the  boundary  layer,  so  that  results  were  inaccurate.  More¬ 
over,  in  figure  23  for  inviscid  flow,  the  results  of  equation  110  were 
less  accurate  than  those  of  equation  115.  From  the  analysis  of  runs  1 
through  13: 

1.  For  inviscid  flow,  the  differencing  technique  of  section  4.3.2  and 
stabilizing  term  of  equation  115  was  accurate. 

2.  For  viscous  flow  and  the  differencing  technique  of  section  4,3.2, 
the  stabilizing  term  of  equation  116  yielded  accurate  results  in  the 
boundary  layer. 

3.  Other  numerical  techniques  and  coding  used  in  conjunction  with 
items  1  and  2  were  satisfactory. 

4.  Solutions  involving  equation  110  with  Cj  =  =  3  converged  rapidly 

and  eliminated  perturbations  but  were  inaccurate  in  the  boundary  layer, 

5.  Nodal  spacing  was  critical  for  stability  and  ac-iuracy. 

By  incorporating  these  ideas  into  the  specifications  of  runs  14  and 
15,  a  satisfactory  solution  was  obtained  for  viacous  flow  around  a 
hemisphere-cylinder.  The  small  Reynolds  number  of  4000  was  selected 
to  reduce  the  number  of  cycles  required  for  convergence.  Experimen¬ 
tation  revealed  that  stability  was  enhanced  by  increasing  the  number  of 
nodes.  The  number  of  nodes  was  384  in  runs  12  and  13  and  was  841  in 
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runs  14  and  15.  In  run  14  with  the  stabilizing  term  of  equation  110  with 
Cj  =  C^  =  3,  the  perturbations  of  initial  values  were  not  apparent  in  the 
results  at  K  =  1600.  Using  the  output  tape  of  run  14  for  initial  values 
in  run  15,  the  results  converged  within  800  cycles.  As  illustrated  by 
figure  29,  the  effect  of  the  stabilizing  term  was  small,  even  in  the 
boundary  layer.  Results  of  run  15  are  presented  in  section  5.5.3. 

In  applying  equation  llo  to  other  viscous  systems,  acceptable 
values  of  Cg,  Cg,  and  may  be  determined  by  numerical  experi¬ 
ments  in  which  relative  magnitudes  of  and  are  compared.  Addi¬ 
tional  numerical  experiments  may  be  desirable  to  optimize  the  number 
of  nodes  and  nodal  spacing. 
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Figure  23*  Computed  v-velocUtee.  Near  Junction  at  1  -  16, 
Inviecid  Flow  around  Hemlephere-Cyllnder.  M  »  4* 
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Vlcous  Term  tp. 


Figure  25,  Viscous  Term  <p-  and  Stabilizing  Term  at  I  =  8 

and  K  =  2000.  Run  ll.  Hemisphere  at  M  =4  and  Re  =  10  . 

o  o 

Computed  by  Section  4,  3,  2  with  Equation  116, 
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viscous  Term  © . 

4 

Figure  26.  Viscous  Term  ©  and  Stabilizing  Term  at  I  -  8 

and  K  =  2000.  Run  ll.  Hemisphere  at  M  =4  and  Re  =  10  . 

o  o 

Computed  by  Section  4.  3. 2  with  Equation  116. 
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5.4  Accuracy 

5.4. 1  Comparisons  with  Measurements.  The  accuracy  of  results 

could  not  be  determined  analytically  and  was  estimated  ly  comparisons 

with  computed  and  measured  results  of  others.  Measurements  of  bow 

wave  coordinates  and  body  surface  pressures  are  available  in  [4]  at 

5  i 

M^  =  2  and  4  and  =  4.  93  x  10  and  1.28  x  10  respectively.  Steady 
solutions  of  viscous  flows  for  such  high  Reynolds  numbers  would  require 
excessive  computing  time  and  were  not  attempted.  It  was  assumed  that 
computed  inviscid  results  of  runs  6  and  16  should  agree  with  these 
measurements.  Measured  and  computed  values  were  not  at  the  same 
points  and  differences  of  results  were  estimated  from  graphs.  Computed 
and  measured  results  are  shown  In  figures  20  and  21  on  pages  96  and  97 
for  M^  =  4  and  in  figures  30  and  31  on  pages  108  and  109  for  M^  -  2. 

Computed  and  measured  [75]  u/u^<*velocitles  in  the  boundary  layer 
at  several  locations  are  displayed  in  figures  32  through  34.  The  agree¬ 
ment  verifies  the  accuracy  of  the  numerical  techniques  for  the  boundary 

layer.  The  central  angle  9  of  computed  and  measured  results  is 

c 

slightly  different  because  of  nodal  requirements. 

Tabulated  average  and  maximum  differences  s  e  presented  in 
Table  3.  In  these  comparisons,  maximum  differences  as  high  as  13. 7 
percent  were  noted  at  selected  points  where  either  measurements  or 
calculations  were  difficult  to  obtain.  However,  the  average  difference 
of  measured  and  computed  results  was  5.4  percent  or  less.  Considering 
the  probable  error  of  the  measurements,  agreement  of  results  is  very  good. 
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Table  3.  Percentage  Differences  of  Computations  and  Measurements 


Variable 

Max.  % 
Diff. 

Avg.  % 
Diff. 

Location 
Max.  Diff. 

M 

o 

Re 

o 

Run 

Experi¬ 

ment 

Figure 

Pb 

11.7 

1.9 

X=1.57 

4 

mv 

6 

[4] 

22 

Pb 

1.0 

.4 

X=1.02 

2 

mv 

16 

’C4] 

31 

u> 

w 

9.9 

4.5 

X=1.91 

4 

INV 

6 

[4] 

23 

03 

w 

6.2 

4.8 

X=1.02 

2 

INV 

16 

[4] 

30 

d 

w 

4.9 

x=o 

4 

INV 

6 

[4] 

23 

d 

w 

5.3 

x=o 

2 

INV 

16 

[4] 

30 

u/u 

e 

4.4 

2.1 

0)=.  0144 

4 

10^ 

17 

L'75] 

32 

u/u 

e 

13.7 

5.4 

0)=.  006 

4 

10^ 

17 

[75] 

33 

u/u 

e 

5.7 

2.5 

Ci)=.  016 

4 

10^ 

17 

[75] 

34 

Wave  Radius 


— •  Run  17,  Section  4.  3, 2,  Equation  116, 

e  =32.8®,  u  =.8368 
c  e 

©  Experl3nentC75],  9^  =  30° 


Figure  32.  Boundary  Layer  Velocity  Ratios  around  Hemisphere 
at  M  =2  and  Re  =  10^. 


Figure  33.  Boundary  Layer  Velocity  Ratios  around 

Hemisphere  at  M  ~  2  and  Re  =  10^. 

o  o 
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Figure  34.  Boundary  Layer  Velocity  Ratios  around 
Hemiei^ere  at  =  2  and  Re^  =  10  . 
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5.4.2  Comparison  with  Computed  Resuits  of  Others.  Prior  to  this 
investigation,  computed  results  were  not  available  for  viscous  flow 
around  a  hemisphere  or  hemisphere-cylinder.  Therefore,  comparisons 
of  computed  results  were  restricted  to  inviscid  flow  systems.  Compar¬ 
isons  were  made  with  numerical  results  of  Aungler  [1]  and  analytical 
results  of  Belotserkovskii  [5]  at  =  2  and  4.  The  latter  analytical 
results  were  obtained  by  a  method  of  integral  relations  that  is  discussed 
in  section  2.  5.  A  summary  of  differences  in  results  is  presented  in 
table  4.  Entropy  along  the  body  surface  was  constant  for  run  6  and 
reference  1.  The  difference  of  those  entropies  was  0.2  percent.  Body 
surface  entropy  values  were  not  presented  in  [1]  for  =  2.  For 
inviscid  flow,  body  surface  entropies  were  calculated  theoretically. 
Comparisons  with  those  values  are  also  included  in  table  4.  As  indi¬ 
cated  by  tables  3  and  4,  the  results  of  runs  6,  16  and  17  were  in  agree¬ 
ment  with  those  of  [1],  [4],  [5]  and  [75].  These  solutions  demon¬ 
strated  that  the  computer  code  developed  in  this  investigation  is 
satisfactory  for  inviscid  and  viscous  compressible  flows  around  axi- 
symmetrical  bodies. 
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Table  4.  Percentage  Differences  of  Computed  Results  of  this 
Investigation  and  Others. 


Variable 

Max.  % 
Diff. 

Avg.  % 
Diff. 

Location 
Max.  Diff. 

Run 

Reference 
or  Theory 

Pb 

wO 

«0 

2 

16 

cn 

Pbi 

s=0 

»0 

4 

6 

[1] 

CO 

w 

4.8 

3.2 

X=1.28 

2 

16 

[1] 

CO 

w 

3.  1 

1.8 

X=  .84 

4 

6 

cn 

d 

w 

3.4 

x=o 

2 

16 

[1] 

d 

w 

2.2 

x=o 

4 

6 

[1] 

d 

w 

.  0.8 

x=o 

4 

6 

[53 

% 

0.2 

0.2 

x=o 

4 

6 

[1] 

^b 

1.4 

2 

16 

Theory 

«b 

0.4 

4 

6 

Theory 

5.  5  Results 

5.5. 1  Introduction.  This  section  includes  solutions  of  inviscid  and 
viscous  flows  around  hemispheres  and  hemisphere -cylinders.  Prior 
to  this  investigation,  time  dependent  solutions  of  viscous  compressible 
flow  about  hemispheres  and  hemisphere-cylinders  were  not  available 
in  the  open  literature. 

5.5.2  Results  for  a  Hemisphere.  Specifications  of  runs  6,  11,  18, 

19  and  20  for  hemispheres  are  presented  in  table  2.  Each  run  was 
preceded  by  numerical  experiments  to  determine  satisfactory  param¬ 
eters  for  the  stabilizing  terms  of  equations  115  and  116.  As  in 
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section  5.3.3,  comparisons  vtere  made  of  stabilizing  and  viscous  terms. 

Runs  6  and  11  at  M  =4  with  inviscid  and  viscous  flows  were  a  part  of 

o 

the  development  phase.  Their  results  are  described  in  section  5.  3.  3. 

5 

Run  20  for  M  =4  and  Re  =10  was  difficult  because  of  the  thin  bound- 
o  o 

ary  layer.  Since  the  total  number  of  nodes  was  large  and  Ar  was  small, 
the  computing  time  for  run  20  was  greater  than  in  other  runs.  Results 
were  nearly  steady  at  K  =  3500. 

For  =  4,  the  standoff  distance  shown  in  figure  35,  page  117, 

decreases  as  Reynolds  number  increases.  In  figure  36,  there  is  a 

decrease  in  body  surface  pressure  for  an  increase  in  Reynolds  number. 

5 

In  that  figure,  inviscid  flow  and  Re^  =  10  are  represented  by  a  single 
line  for  X  <  .8.  The  boundary  layer  displacement  thicknesses  of 
figure  37  were  computed  from 


(151) 


where:  m  =  pu 

e  refers  to  the  edge  of  the  boundary  layer. 

The  limit  was  determined  from  a  graph,  such  as  figure  38,  of  u- 

velocity  versus  co  at  the  point  where  the  boundary  layer  effects  seemed 

4  5  * 

negligible.  For  Re^  =  10  and  10  ,  there  are  oscillations  in  6  near 

the  stagnation  point  where  the  boundary  layer  is  thin  because  of  small 

oscillations  in  u-velocity  near  the  edge  of  the  boundary  layer.  For  the 

3  3 

relatively  thick  boundary  layers  of  =  10  and  4x10,  these 
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oscillations  had  a  minor  effect  on  d  .  In  figures  38  through  40,  u- 
velocities,  pressures  and  temperatures  are  presented.  The  decrease 
in  boundary  layer  thickness  with  increasing  Reynolds  numbers  is 
apparent. 
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Body  Pressure 


Figure  36.  Body  Surface  Pressures  of  Hemisphere  at  =  4. 
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5.5,3  Results  for  a  Hemisphere -Cylinder.  Based  on  the  numerical 
results  of  section  5.3.3,  the  specifications  of  runs  14  and  15  of  table  2 
were  selected  for  a  solution  of  viscous  compressible  flow  around  a 
hemisphere-cylinder.  In  run  14  for  K  s  1600,  the  stabilizing  term  of 
equation  110  was  used  to  remove  initial  perturbations.  By  KL  =  1600, 
a  favorable  pressure  gradient  had  been  established  in  the  afterbody 
region,  and  the  bow  wave  had  moved  close  to  its  steady  location.  The 
results  of  run  14  at  KL  =  1600  were  used  as  input  to  run  15  in  which  the 
stabilizing  term  was  that  of  equation  116.  In  run  15  at  KL  =  800,  the 
solution  was  steady  and  seemed  accurate  in  the  boundary  layer  and 
other  parts  of  the  field.  Results  of  run  15  for  KL  =  800  are  presented 
in  figures  41  through  56  at  the  end  of  this  section  on  pages  127  to  142. 
Contours  of  Mach  number,  temperature,  pressure  and  entropy  are 
presented  in  figures  41  through  44.  The  contours  are  consistent  with 
each  other.  Near  the  downstream  boundary  and  bow  wave,  there  are 
oscillations  in  the  results.  The  contours  for  Mach  number  and  temper¬ 
ature  in  figures  41  and  42  are  similar.  Boundary  layer  growth  along 
the  body  is  apparent.  Mach  number  contours  outside  of  the  boundary 
layer  in  the  forebody  region  are  similar  to  those  of  inviscid  flow  in 
[52].  In  figure  43,  as  expected,  there  is  an  expansion  fan  near  the 
junction  of  the  hemisphere  and  cylinder.  Entropy  contours  of  figure  44 
were  smooth  and  resembled  streamlines  outside  of  the  boundary  layer. 
Along  the  isothermal  body  surface  Tj^  =  4.  2  and  entropy  increases 
because  pressure  decreases. 
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Boundary  layer  displacement  and  momentum  thicknesses  of  figure 
45  seem  reasonable.  There  is  a  substantial  increase  of  boundary  layer 
thickness  along  the  afterbody.  Figures  46  through  49  present  u- 
velocities,  v-velocities  and  temperatures  versus  o)  at  several  I  indices. 
The  coordinate  a)  was  chosen  instead  of  (o/o)^  to  display  growth  of  the 
boundary  layer.  Variations  in  u-velocities  and  temperatures  are 
smooth.  At  I  =  16  near  the  junction  in  figure  46,  there  is  a  small  de¬ 
crease  in  u-velocity  near  the  edge  of  the  boundary  layer.  Near  the 
forebody,  variations  in  v-velocities  of  figure  47  are  smooth.  As  shown 
in  figure  48,  they  uscUlate  in  the  afterbody  region  where  they  are  small 
and  sensitive  to  perturbations.  Although  undesirable,  those  oscillations 
were  not  considered  detrimental  to  the  accuracy  of  other  variables.  In 
figure  49i  temperature  gradients  <  0  and  their  absolute  values 

increase  with  increasing  I  indices.  Since  body  temperature  is  greater 

than  adiabatic  stagnation  temperature,  heat  transfer  is  from  the  body 
9  T\ 

and  •r“|  should  be  negative. 

aw/b 

Specifications  of  runs  15  and  19  in  table  Z  for  a  hemisphere- 
cylinder  and  hemisphere  were  identical.  As  indicated  by  table  5, 
results  of  both  runs  in  the  forebody  region  are  about  the  same. 


Table  5.  Percentage  Differences  of  Computed  Results  of  Hemisphere 
and  Hemisphere -Cylinder.  Runs  15  and  19,  M  =  4,  Re  =  4000. 


Max.  % 


Av^.  % 


Location  of  Maximum 


Difference 

Difference 

1  Differenc 

3.3 

.6 

o 

II 

1.3 

.  1 

X  =  1.28 

3.  0 

.3 

X  =  1.28 

1.2 

.4 

X  =  1.28 

.4 

o 

II 

A  conclusion  was  that  satisfactory  results  for  the  forebody  regio;' 
hemisphere -cylinder  may  be  obtained  by  considering  the  hemisphere 
alone.  The  conclusion  may  be  true  for  other  shapes  and  conditions,  but 
additional  results  are  needed  to  strengthen  this  conclusion. 

In  run  21,  the  stabilizing  term  was  zero  and  results  of  run  15  at 
KL  =  800  were  used  for  initial  values.  In  run  21,  results  were  steady 
at  KL  =  1600.  Typical  results  of  runs  15  and  21  for  the  forebody  region 
are  shown  in  figures  50  through  52.  The  maximum  difference  is  9  per¬ 
cent.  Average  differences  of  all  variables  is  less  than  2  percent.  The 
variation  in  pressure  versus  K  index  in  the  afterbody  region  for  various 
stabilizing  terms  is  shown  in  figure  53  at  I  =  19  and  J  =  26  where  differ¬ 
ences  of  runs  15  and  21  were  a  maximum.  The  large  differences  at  that 
point  were  due  to  the  spatial  oscillations  in  the  results  of  run  21.  u- 
velocities,  v-veloclties  and  pressure  for  run  15  and  21  at  I  =  19  are 
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shown  in  figures  54  through  56.  Oscillsitions  of  the  steady  results  of 
run  21  are  a  maximum  near  the  bow  wave  and  are  insignificant  in  the 
boundary  layer.  Where  there  are  no  significant  oscillations  of  results, 
the  differences  between  runs  15  and  21  are  less  than  3  percent.  The 
oscillations  of  run  21  without  stabilizing  terms  were  about  the  smooth 
values  of  run  15.  The  agreement  of  the  results  of  runs  15  and  21  indi¬ 
cates  that  the  accuracy  of  the  results  of  run  15  were  not  seriously 
affected  by  the  addition  of  stabilizing  terms. 
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Figure  41.  Mach  Number  Contours  for  Hemisphere -Cylinder  at  M  =4  anc 
Run  15.  Computed  by  Section  4.  3.2  and  Equation  116. 


Figure  4Z,  Temperature  Contours  for  Hemisphere-Cylinder  at  =  4  and  ~  4000. 
Run  15.  Computed  by  Section  4.  3.2  and  Equation  116. 


v-velocity 

Figure  47.  v-velocitles  in  Forebody  Region  of  Hemisphere-Cylinder 
at  M  =4  and  Re  =  4000,  Run  15.  Computed  by  Section  4.  3.  2  and 
Equation  U6.  ° 


;7J/  ^7^K7;l  .  --»v^ 


—  Run  ’.5,  Equation  116 
o  Run  21,  =  0 


U  l.U  1  11. 

u-veloclty 

Figure  50.  u-velocltles  In  Forebody  Region  of  Hemisphere -Cylinder 

at  I  =  8  with  and  without  Stabilizing  Terms.  Computed  by  Section 

4.  3.  2,  M  =4  and  Re  =  4000. 
o  o 
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v-veloclty 

Figure  51.  v-velocitles  In  Forebody  Region  of  Hemisphere-Cylinder 

at  I  =  8  with  and  without  Stabilizing  Terms.  Computed  by  Section 

4.3.2.  M  =  4  and  Re  =  4000. 
o  o 
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Run  15,  Equation  116 
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Figure  56.  Pressures  Ir  Afterbody  Region  of  Hemisphere-Cylinder 

at  I  =  19  with  and  without  Stabilizing  Terms.  Computed  by  Section 

4.3.2.  M  =  4  and  Re  =  4000, 
o  o 
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5,  6  Summary  of  Results 

Considerable  analyses  and  computations  preceded  the  successful 
solution  of  viscous  flow  around  a  hemisphere -cylinder  to  develop  tech¬ 
niques.  Those  analyses  and  computations  seemed  to  indicate: 

1.  The  body  related  coordinates  of  section  3.3.2  are  applicable  to 
most  blunt  bodies. 

2.  The  nonlinear  coordinate  transformations  of  section  3,3.3  permit 
adequate  resolution  in  the  boundary  layer  without  an  excessive  total 
number  of  nodes. 

3.  To  obtain  st'*ady  results,  initial  values  should  be  as  accurate  as 
possible  to  enhance  stability  and  reduce  computer  time. 

4.  Techniques  for  boundary  values  at  the  body  surface  and  downstream 
boundary  have  a  major  effect  on  stability  and  ac  curacy. 

5.  Precise  application  of  the  method  of  characteristics  was  the  only 
technique  that  was  satisfactory  for  wave  fitting. 

6.  The  differencing  technique  of  section  4.  3.  2  with  stabilizing  terms 
should  be  applicable  to  a  wide  variety  of  axisymmetrical  invlscid  or 
viscous  flows. 

The  most  valuable  result  of  this  investigation  is  the  complete  set 
of  numerical  techniques  of  section  4  for  axisymmetrical,  viscous, 
compressible  flows  around  blunt  bodies.  Results  are  presented  for 
viscous,  compressible  flow  around  a  hemisphere- cylinder  at  free 
stream  Mach  and  Reynolds  numbers  of  4  and  4000  and  for  flows  around 


hemispheres  at  Mach  numbers  of  2  and  4  and  Reynolds  numbers  of 
10^,  4x10^,  10^  and  10^. 
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6.  CONCLUSIONS 

1.  A  complete  set  of  numerical  techniques  for  axisymmetrical , 
viscous,  compressible  flows  around  blunt  bodies  was  developed. 

2.  An  accurate  solution  was  obtained  for  viscous  compressible  flow 
around  a  hemisphere-cylinder  for  freestream  Mach  and  Reynolds 
numbers  of  4  and  4000  respectively. 

3.  Numerical  stability  was  achieved  with  stabilizing  terms,  but  their 
necessity  was  not  established. 

4.  Approximate  results  can  be  obtained  for  the  forebody  region  of  a 
hemisphere-cylinder  with  a  free  stream  Mach  number  of  4  and 
Reynolds  number  of  4000  by  solving  the  corresponding  system  of  a 
hemisphere  alone. 
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